Fortran 90
for the Fortran 77 programmer
Copyright C 1993, Bo Einarsson and Yurij Shokin
Permission is granted to copy and/or print this file as long
as the copyright notice and this permission is included on
all copies. This file is also available as one book in Swedish
and one in Russian. A complete textbook in Swedish on Fortran 90
has been prepared, and is available both electronically and in
printed form.
The pages are at present numbered as in the Swedish version of
24 January 1994. This regrettably implies that some pages are
longer than a standard page on a printer.
Authors's addresses
Bo Einarsson Yurij I Shokin
National Supercomputer Centre Institute of Computational Technologies
Prospekt Lavrentyeva 6
University of Linkoping Russian Academy of Sciences
Siberian Division
S-581 83 LINKOPING SU-630090 NOVOSIBIRSK 90
SWEDEN RUSSIA
Tel. + 46 13 281432 + 7 383 235 00 50
Fax. + 46 13 282535 + 7 383 235 12 42
Email: boein@nsc.liu.se shokin@ict.nsk.su
WWW: http://math.liu.se/~boein/
This version is dated 23 January 1995.
--------------------------------------------------------------------
Dieses Skriptum ist so formatiert, daß es als HTML-File gleichzeitig mit
lpr -Pf68umps
als Postscript-File mit korrektem Seitenumbruch ausgedruckt werden kann,
z.B. mit Netscape Print... mit Paper Size A4. Die Schriftgröße muß mit
Options General_Preferences Fonts Size:12.0
eingestellt sein (Proportional- und Fixed Font)
--------------------------------------------------------------------
Preface v
Introduction 1
1. Notables in the transition from Fortran 77 to Fortran 90 3
2. Specifications 4
3. The layout of a program (free form and fix form) 6
4. Format 8
5. The same source code for Fortran 77 and Fortran 90 ? 9
6. Control statements 9
7. Program units 11
8. Keyword arguments and default arguments 14
9. Recursion 16
10. Generic routines 17
11. The use of arrays and array sections 19
12. Pointers 21
13. The new precision concepts 24
14. Additional problems at the transition 24
15. Use of program libraries 26
16. Peculiarities in the language Fortran 90 27
Answers and comments to the user exercises 27
References 38
Appendices
1. Fortran and Pascal 41
2. Summary of Fortran 77 statements 43
3. Summary of the new features in Fortran 90 49
4. Backward and Forward compatibility 57
5. Intrinsic functions and subroutines in Fortran 90 58
6. Fortran 90 from NAG 72
7. Other Fortran 90 systems 75
8. Fortran 90 properties in Cray CF 77 76
9. The historical development of Fortran 77
10. High Performance Fortran 83
11. Explanations of certain terms 93
12. Argument association 97
13. Index 100
Marketing information, the book 110
Marketing information, the authors 111
----------------------------------------------------------------------------
- v -
Preface
This book is written in order to ease the transition from the very
common and popular programming language Fortran 77 to the more modern
Fortran 90. This transition uses the fact that Fortran 77 is a pure
subset of Fortran 90. There are, however, two very important reasons to
go over to as much Fortran 90 as possible. One is that it includes new
and powerful constructs, the other is that Fortran 90 gives us more
facilities for correctness checking of the program. This means that
more reliable programs are obtained. During the end of 1994 Fortran 90
will be offered by most computer manufacturers and the language will be
a success in 1995.
It is required that the reader is knowledgeable in Fortran 77. Those
parts of Fortran 90 who already are parts of Fortran 77 are not treated
systematically in full, they are assumed known by the reader. Those who
don't know Fortran 77, can read a tutorial book on Fortran 77, or a
complete textbook on Fortran 90 for example "Programmer's Guide To
Fortran 90" by Brainerd, Goldberg and Adams (McGraw - Hill 1990). Note
especially that Fortran 90 is much larger than Fortran 77 in all
respects. Therefore it is difficult to describe it as short as we do it
here. All examples in this book have been run on a Sun SPARC (and then
also on the DEC station Ultrix and the IBM PC). All tests have been
done using the Fortran 90 system from NAG. It is recommended that the
reader has access to a Fortran 90 system.
The examples in the book are considered to illustrate the programming
technology and the statements being discussed.
The purpose is, however, not to give an optimized application
program. This is especially true for the sections supplying comments to
the exercises. Please note that some of the later examples are very
complete with respect to interfaces and specifications that are required
at the use of functions and subroutines.
The proposed extension HPF or High Performance Fortran, primarily
intended for parallel computers, is described in an appendix.
We wish to thank John Reid for his assistance with the interpretation
of some of the new concepts in Fortran 90.
Views, suggestions and corrections are most welcome for the next
edition of this tutorial.
[ Translation remarks and queries are indicated in this way. ]
Linkoping and Novosibirsk, August 1993
Bo Einarsson and Yurij Shokin
----------------------------------------------------------------------------
- 1 -
Introduction
The programming language Fortran is the dominating language for
technical and scientific applications. It was developed originally at
IBM in 1954 under the supervision of John Backus. Already in 1958 the
language was expanded to Fortran II, which included subroutines,
functions and common blocks. Other manufacturers came later with
Fortran compilers for their machines and in 1962 IBM introduced the
extended Fortran IV, which was the base for the 1966 agreed American
standard ANSI X3.9 - 1966. The International Standardization Committee
agreed on three different levels of Fortran in the document ISO 1539 -
1972.
In April 1978 the American Standardization Committee ANSI approved a
new standard for Fortran, it is well known with the name Fortran 77 in
order to distinguish the new standard ANSI X3.9 - 1978 from the old ANSI
X3.9 - 1966, which sometimes nowadays is called Fortran 66. Fortran 77
was very welcome in that it increased the capabilities of writing
Fortran programs in a structured way, and that it standardized many of
the extensions that different manufacturers had put into their
implementations of Fortran.
The International Standard Committee accepted Fortran 77 as a
standard with the document ISO 1539 - 1980. Also the United States
government has accepted Fortran 77 as a Federal standard. The American
Standardization Committee has now decided that Fortran 77 will be a
permanent standard with the new Fortran 90 considered as a different
language within the United States.
Fortran has now been revised considerably. The working group X3J3 at
this time also had some European members, which have had a great
influence. Regrettably there is a naming conflict: internationally the
new work is called Fortran while in America the old version is called
Fortran and new version Fortran 90. In order to avoid confusion it's
therefore advisable to use the name Fortran 77 for the old version and
Fortran 90 for the new version.
The International Standardization Committee has accepted Fortran 90
as a new standard ISO 1539 - 1991. Internationally therefore Fortran 77
is no longer a standard. In September 1992 the American Standardization
Committee accepted Fortran 90 as ANSI X3.198 - 1992.
The purpose of the work on a new standard version was to make Fortran
usable and efficient language for scientific and technical computations
during the nineties. At the moment there is regrettably only a very
small number of compilers available for Fortran 90, but that situation
will soon improve. The new version of the language contains new
facilities for vector and matrix operations, which is very good on
vector processors and parallel processors and also several new methods
to specify precision (not only the REAL and DOUBLE PRECISION),
possibilities to ask for environmental parameters, intrinsic functions
to manipulate floating-point numbers (find an exponent, find the base or
find the mantissa, and to put those parts together into a floating point
number), intrinsic procedures and new specifications for storage and
interfaces. In addition, there is a new improved form for the source
code and more control statements, recursion and dynamically allocatable
arrays are introduced.
----------------------------------------------------------------------
- 2 -
Nothing was removed, therefore Fortran 90 includes both modern and
not so modern methods for essentially the same tasks. This means that
the language is large. The Committee has therefore chosen to call
certain concepts in Fortran 77 obsolete, and to perhaps remove them at
the next revision, see Appendix 4 "BACKWARD COMPATIBILITY".
Note that at the standardization 1977 - 78 some new features were
included like IF...THEN...ELSE...END IF, floating point variables as
index in DO-loop (silly that these were introduced) and also the new
basic data type CHARACTER for text. The extended DO-loop was removed
(which was a very good thing) and the Hollerith constant was removed
(except in FORMAT). That means that there are a few programs that obey
the Fortran 66 standard but do not obey the Fortran 77 standard. Those
things that were removed from the standard in 1978 are included in some
of the implementations of Fortran 90. Extensions are permitted, if they
can be signaled by the system as differences from the standard.
A very important requirement at the introduction of a new or revised
programming language is nowadays that it should permit efficient
compilation and execution, not only on conventional computers but
preferably also on parallel and massively parallel systems. Somewhat
simplified it can be said that Fortran 90 is efficient up to and
including vector processors, but less efficient on parallel systems.
The reason is simply that Fortran 90 does not contain any statements for
division into parallel execution, except what is included automatically
in the array concept.
The proposed extension HPF or High Performance Fortran, primarily
intended for parallel computers, is therefore described in Appendix 10.
This language has the ability to control the allocation of arrays on
parallel systems and include facilities for
Parallel execution
Good performance on any type of computer
Allocating data on the available processors
The main part of this tutorial is sections 1 to 16 on pages 3 to 27
and has the purpose of giving the reader an introduction to Fortran 90.
The understanding is greatly improved if the reader does the included
Exercises. Most of these have complete solutions on pages 28 to 37.
The Appendices 2 and 3 summarize Fortran 77 and the new facilities in
Fortran 90, respectively. These two appendices can together with
Appendix 5 with the intrinsic functions and subroutines be used as a
reference. The Word explanations and the Index can be used for finding
additional information. Those that plan to get a compiler are
recommended to read Appendices 6 and 7.
----------------------------------------------------------------------
- 3 -
Fortran 90
I assume that the reader can write programs in Fortran 77 and wishes
to learn to use new facilities in Fortran 90. All statements in Fortran
77 are explained in Appendix 2 and the summary of the news in Fortran 90
are in Appendix 3.
The greater power in Fortran 90 means that the statements in many
cases have a combined effect, and it is therefore not so useful to only
describe the language statement for statement.
1. Notables in the transition from Fortran 77 to Fortran 90
============================================================
o Everything that is in Fortran 77 is also in Fortran 90.
o Many new statements have been added, some replacing older
statements.
o Many new statements have been added and give new possibilities.
o There are now two forms of the source code. The old source code
form, which is based on the punched card, and now called fixed
form and the new free form.
These two forms may not be mixed independently but must be
separated at the compilation; an old program, (one or several
program units in the form of a main program, subroutines or
functions), which is compiled in free form, has the possibility to
give a different result compared with earlier when it was compiled
with fixed form, compilation errors are possible.
In the same program we can mix program units written in fixed
form and free form, but each unit must be only in one form and
at the compilation both forms may not usually be in the same
source file. Some systems using a special directive may permit
program units in different form in the same file. The directive
then tells the compiler which form is valid.
o Some old statements are to be avoided.
o It is possible to mix old and new statements, but it is advised
to try to be consistent, which means to use either the old
or the new form for the statement.
o Note that when you switch from one compiler to another new
errors may occur, because the old compiler was not as strict as
the new one. Normally, a new compiler discovers some old errors
that were not found earlier. In the specification for Fortran 90
it is required that errors are to be found already at compilation
if this is at all possible.
Exercises
---------
(1.1) Compile and run some of your own small programs with
Fortran 90 using fixed form. Get some instructions in
Appendix 6, NAG's Fortran 90.
(1.2) Modify the program you have just compiled with exchanging
comments starting with C or * for comments starting with !
and try to use Fortran 90 free form.
----------------------------------------------------------------------
- 4 -
(1.3) What happens if the following small program is run using
fixed form or using free form?
LOGICAL L
L = .FALSE.
IF (L) THEN Z = 1.0
ELSE Y = Z ENDIF
END
2. Specifications
=================
My favorite statement in Fortran 90 is IMPLICIT NONE, which means
that the implicit declaration is no longer used. This statement has
been available in some implementations of Fortran 77. The use of that
statement means that the probability of errors because of wrongly spelt
variable names is drastically reduced.
The first difference regarding specification of variables is that
these can now be put together in one statement for each variable. Using
Fortran 77 you can declare, for example, as follows
REAL A, B, C
PARAMETER (A = 3.141592654) DIMENSION B(3) DATA B /1.0, 2.0, 3.0 /
DIMENSION C(100) DATA C /100*0.0/
that is one variable can occur in several lines. Using Fortran 90
you can instead write
REAL, PARAMETER :: A = 3.141592654
REAL, DIMENSION(1:3) :: B = (/1.0, 2.0, 3.0 /)
REAL, DIMENSION(1:100) :: C
DATA C /100*0.0/
The difference is not so large here but it's much larger in more
complicated examples with many variables, especially since in Fortran 90
you have access to more properties. It is astonishing that you can not
specify that all the values of the vector C should be zero directly in
the specification. The Fortran 90 has no support for the repetition
factors in a statements like
REAL, DIMENSION(1:100) :: C=(/100*0.0/).
Neither can you write in the following manner
C = (/0*I, I=1,100/)
If there are variables of different types there is an intrinsic
function that shows what the exact subtype of the variable used. This
function is called KIND (x). This function can also be used for
particular types of subtypes or kinds:
KIND (0) integer
KIND (0.0) floating point number
KIND (.FALSE.) logical variable
KIND ("A") string of characters
There is an intrinsic function SELECTED_REAL_KIND, which returns the
kind of the type REAL that has a representation with (at least) the
precision and the exponential range requested. For example,
-------------------------------------------------------------------------
- 5 -
the function SELECTED_REAL_KIND (8,70) is that kind of REAL that has
at least 8 decimal digits accuracy and permits an exponent between
10**-70 and 10**70. The corresponding function for integers is called
SELECTED_INT_KIND and of course has only one argument. With the
choice SELECTED_INT_KIND (5) all integers between (but not including
the limits) -100 000 and +100 000 are permitted. The kind of a type
can be given a name.
INTEGER, PARAMETER :: K5 = SELECTED_INT_KIND(5)
This kind of integers can be used in constants according to the
following line
-12345_K5
+1_K5
2_K5
which is a rather unnatural specification, after the value we have to
give an underscore _ followed by the name of the kind.
Use of variables of the new integer type can be declared in a nicer
way
INTEGER (KIND=K5) :: IVAR
The corresponding is true for floating-point variables, if we first
introduce a high-precision kind LONG with
INTEGER, PARAMETER :: LONG = SELECTED_REAL_KIND(15,99)
then we get the floating-point kind with at least 15 decimal digits
accuracy and with an exponent range from 10**-99 to 10**+99. The
corresponding constants are obtained as
2.6_LONG
12.34567890123456E30_LONG
and variables are declared with
REAL (KIND=LONG) :: LASSE
The old type conversions INT, REAL and CMPLX have been extended with
the functions
INT(X, KIND = K5)
which converts a floating-point number X to an integer of the
kind K5, if Z is a complex number with
REAL(Z, KIND(Z))
you get it converted to a floating-point number of the real type
and of the same kind of Z (that is of course the real part of Z).
Double precision is not included in the new Fortran 90 in any other
way than in the "old" Fortran 77, but it is assumed that the compiler
supports the double or quadruple precision that may be available in the
hardware. You can then define a suitable kind of the REAL, named DP or
QP. You can of course use the old concept of DOUBLE PRECISION.
----------------------------------------------------------------------
- 6 -
The reason for this rather cumbersome way is that it is not desirable
to have too many compulsory precisions (for example single, double,
quadruple, perhaps for both the cases REAL and COMPLEX) and also that
the old concept DOUBLE PRECISION did not give a specified machine
accuracy. Now you can relatively easily specify both which precision
and which range of exponent you wish to use. Additional information
about the kind is given in Appendix 6, where the different data types
and their normal kinds on both DEC and SUN and also on the IBM PC is
given in the NAG system for Fortran 90.
Exercises
---------
(2.1) What does the specification LOGIC ALL mean?
(2.2) Specify a constant K with the value 0.75.
(2.3) Specify an integer matrix PELLE with 3 rows and 4 columns.
(2.4) Specify a floating-point number which corresponds to the
double precision on an IBM and a single precision on Cray.
(2.5) Specify some variables of the type above.
(2.6) Specify some constants of the type above.
(2.7) Is the following specification correct?
REAL DIMENSION(1:3,2:3) :: AA
(2.8) Is the following specification correct?
REAL REAL
(2.9) Is the following specification correct?
COMMON :: A
3. The layout of a program (free form and fix form)
===================================================
Sometimes we require more than one line for a statement
Print, * 'This is a long output line this is the first part',&
' and this is the second part!'
Nowadays, in the free form, we continue a line with the symbol "&"
(called ampersand), i.e. with the sign & at the end of the old line
instead of an almost arbitrary character in column 6 of the new line.
With the compiler we are now using it is possible to include the Swedish
characters in character strings and in comments.
[ PLEASE ADD ABOUT RUSSIAN CHARACTERS IN THE RUSSIAN VERSION /Bo]
Sometimes a certain identifier or a certain numerical number has not
sufficient space on a line. We can then interrupt the identifier
anywhere with the character "&" and then on the next line give a new "&"
as the first non-blank character. You continue then directly from this
"&" without any extra blank. The character "&" therefore works as a
kind of delimiting or syllabification sign. You can write
PI = 3.141592653589793
----------------------------------------------------------------------
- 7 -
or you can write equivalently
PI = 3.14159265&
&3589793
Please note that comment lines can not be continued. The reason for
this is that also in a comment line the sign "&" is treated as belonging
to the comment. However, you can add comment lines also inside the
continuation lines. Note that it is the final "&" of a line that
informs of a continuation line, it is therefore possible to write a text
string of characters including the character "&".
Sometimes you may wish to do it in the opposite way, to have several
statements on the same line. This is being done with the use of
semicolon and with the exclamation mark we can include a comment on the
line.
A = 0.0 ; B = 1.0 ; C = 2.0 ! Initialization
A line may include up to 132 characters, a statement up to 2640
characters, blanks included.
Note that on using the free form blanks are significant, as in my
favourite example:
DO 25 I = 1. 25
This gives a compilation error since the compiler does not find a
comma between the lower and upper the limits, but the compressed version
DO25I=1.25
gives the same result as the non-compressed form and as in Fortran 77
or using the old form (fix form) of Fortran 90, namely that the variable
DO25I is being assigned the value 1.25.
Comments are started with "!" (exclamation mark) and ended with the
end of line. The old types of comments introduced with C or * in column
1 are no longer permitted if you use free form, but of course if you use
the fixed form.
Upper case and lower case characters are equivalent except in
character strings.
The above applies to the new free form. Also in the old, column
-oriented or fixed form, we can use a semicolon or exclamation mark
between columns 7 and 72 but you can not continue with "&" or write
comments in columns 1 to 6 or 73 to 80. Exclamation mark in column 1 of
course means a comment line also in the old fixed form. Some
possibilities of longer lines do exist within the old fixed form
(implementation dependent).
Exercises
---------
(3.1) What does the following line mean?
A = 0.0 ; B = 370 ! First variables ; C = 17.0 ; D = 33.0
(3.2) Are the following lines correct according to Fortran 90?
Y = SIN(MAX(X1,X2)) * EXP( -COS(X3)**I ) - TAN(AT&
& AN(X4))
----------------------------------------------------------------------
- 8 -
4. Format
=========
The old Fortran programs have used numbering of the format
statements. However, that doesn't look very well in pure Fortran 90,
where you don't use statement numbers except in a few exceptional cases.
Already in Fortran 77 there was a facility, which however was not used
very much, to use a format variable (of type CHARACTER) instead of a
numbered format. The Format variable was put directly in the
input/output statement.
Now we will show three different ways of doing this assignment. They
all have both good and bad sides. The program follows
PROGRAM FORMAT
IMPLICIT NONE
REAL :: X
CHARACTER (LEN=11) :: FORM1
CHARACTER (LEN=*), PARAMETER :: FORM2 = "( F12.3,A )"
FORM1 = "( F12.3,A )"
X = 12.0
PRINT FORM1, X, ' HELLO '
WRITE (*, FORM2) 2*X, ' HI '
WRITE (*, "(F12.3,A )") 3*X, ' HI HI '
END
In the PRINT statement we use the character string variable FORM1
with the length 11, which is assigned its value in an explicit
assignment statement. The difficulty with this method is essentially
that you have to count manually the number of characters, if it is too
small the NAG compiler will not give a compilation error, but the error
will show up at execution.
In the first WRITE statement we use instead a character string
constant FORM2. The advantage is that with the PARAMETER statement it
is not necessary to give an explicit length of the constant, but it can
be given the length with the statement LEN=*. The bad thing is that we
can not assign the constant a new value.
In the second WRITE statement we use directly an explicit character
string. The difficulty is that the string can not be reused.
Exercises
---------
(4.1) What does this statement give as its output?
WRITE(*, "( HI )")
(4.2) What does the following statement perform?
CHARACTER (LEN=9) :: FILIP
FILIP = '(1PG14.6)'
WRITE(*,FILIP) 0.001, 1.0, 1000000.
----------------------------------------------------------------------
- 9 -
5. The same source code used for Fortran 77 and Fortran 90?
===========================================================
This is possible if you do it in the following way, that is if you
use the new continuation sign "&" at the end of the old line, but in
position 73 so that it doesn't conflict with Fortran 77, and also choose
the "&" sign as the almost arbitrary character in column 6, in order to
get continuation according to Fortran 77. An introductory "&" is "in
principle" neglected by Fortran 90.
program TEST ! column 73
! |
write(*,*) &
& ' test '
end
! |
! column 6
This is really not standard Fortran 77 since neither "&" nor "!" are
in standardized character set. On the other hand, these constructs are
perfect for program segments that are to be included in Fortran 90
program units using the INCLUDE statement (refer to Appendix 3, section
1), sometimes using the old form and sometimes using the new form of the
source code. The program segment is then supposed to be written as if
blanks were significant.
Comments are an incompatibility problem between Fortran 77 and
Fortran 90, but of course not between fixed and free forms of Fortran 90
since the "!" is permitted in both. The exclamation mark is also
permitted in Sun and DEC Fortran 77 (both DEC Station ULTRIX and VAX
VMS) and on the Cray compiler CF77.
Another important condition is of course that the program really
obeys the standard.
6. Control statements
=====================
As conditional or control statements you have IF in many variants
(but essentially not changed from Fortran 77), DO (with some new
variants) and the completely new statement CASE.
The DO-loop should now be ended with the statement END DO and we no
longer need any statement number. In addition, we can use the statement
EXIT to jump out of the DO-loop and CYCLE in order to go to the next
iteration of the present DO-loop. A DO-loop can be assigned a name,
which is done by giving the name before the DO and followed by a colon.
In addition the final END DO can be followed by the name of the DO-loop.
SUMMA = 0.0
ADAM : DO I = 1, 10
X = TAB(I)
EVA : DO J = 1, 20
IF (X > TAB(J)) CYCLE ADAM
X = X + TAB(J)
END DO EVA
SUMMA = SUMMA + X
IF (SUMMA >= 17.0) EXIT ADAM
END DO ADAM
----------------------------------------------------------------------
- 10 -
In the example above, the execution of the inner loop will be
interrupted with a jump to the next cycle of the outer loop, and thus
the variable sum or SUMMA will not be increased, if X is greater than
the given table value. As soon as the sum is at least 17 the outer loop
is also interrupted.
If no name is given in the EXIT or CYCLE statements the present inner
loop is automatically used. With the present inner loop I mean the one
where the EXIT or CYCLE statements being executed are. These statements
then replace the GOTO to the final statement, which was often used in
the old DO-loop. This final statement usually was an CONTINUE
statement.
Also an IF statement can be given a name. In that case the
corresponding END IF ought to be followed by that name.
A new construct in standard Fortran is CASE. It however appeared in
many Fortran dialects before. It can choose a suitable case for a
scalar argument of type INTEGER, LOGICAL or CHARACTER. A simple example
is based on an integer IVAR.
SELECT CASE (IVAR)
CASE (:-1) ! all negative numbers
WRITE (*,*) 'Negative number'
CASE (0) ! zero case
WRITE (*,*) ' Zero'
CASE (1:9) ! one-digit number
WRITE (*,*) ' Digit ', IVAR
CASE (10:99) ! two-digit number
WRITE (*,*) ' Number ', IVAR
CASE DEFAULT ! all remaining cases
WRITE (*,*) ' Too great number'
END SELECT
It is not permitted with overlapping arguments. This means that one
single argument may not satisfy more than one of the cases of CASE. The
default case does not have to be included. If no valid case is found
the execution will continue with the first statement after the END
SELECT. I recommend that you include a DEFAULT and then give an error
message if an argument has a not permitted value.
It is recommended to use the CASE instead of an assigned or computed
GOTO statement, or an arithmetic IF-statement.
Exercises
---------
(6.1) Write a CASE-statement that performs three different calculations
depending on whether the variable is negative, zero or
any of the first odd prime numbers (3, 5, 7, 11, 13) and
performs nothing in all other cases.
(6.2) Write a DO-loop that adds the square roots of 100 given
numbers, but skips negative numbers and concludes the addition
if the present value is zero.
----------------------------------------------------------------------
- 11 -
7. Program units
================
In addition to the four old program units: PROGRAM (that is the main
program), SUBROUTINE, FUNCTION and BLOCK DATA, the new concept MODULE
has been added, as well as some new things in the old units. I repeat
that subprogram is the common concept for both SUBROUTINE and FUNCTION.
I wish to repeat that under Fortran 77 all program units are
essentially on the same level, even if the main program logically is
superior to the subroutines and functions that are called, and even
though you could write a call map that looks like a tree. In reality
the BLOCK DATA is on a higher level and all the other program units are
on the same level, from the Fortran system viewpoint with the main
program just a little above. The exception is the so-called statement
functions with definitions that have to be first in a program unit
directly after the specification and are internal to that unit and
therefore is on a logically lower level. Regrettably, the normal
Fortran programmer does not use statement functions.
The above means that all routine names are on the same logical level,
that means that two different routines, and two different parts of a big
program are not permitted to have the same name. Quite often numerical
and graphical libraries include thousands of functions and subroutines,
and each routine name has to have at most six characters in the name
under old Fortran standards. Therefore, it is very strong risk of a
conflict of names. This problem could be partially solved by the old
statement functions, since these are internal to the respective unit,
and therefore different statement functions can have the same name if
they are in different units. The disadvantage is that they can only
treat what is in only one program line. But they can call each other in
such a way that a later statement subroutine or function can call an
earlier statement function, but of course not the opposite.
Now internal functions and internal subroutines are added, giving a
greater flexibility. They are specified at the end of each program unit
(but not in the BLOCK DATA) after the new command CONTAINS and before
the END. An internal subprogram can have access to the same variables
as the unit it belongs to including the possibility to call its other
internal subprograms. It is written as an ordinary subprogram, but it
is not permitted to have any internal functions or subroutines.
Usual subroutines and functions are as earlier external subroutines
and external functions, but there is now a greater reason for this name
(that is calling them external) than earlier, since now you have also
internal subprograms. Earlier you only had the built in (intrinsic) as
an alternative. In addition, the number of intrinsic functions has
increased very much, and a few intrinsic subroutines have been added.
In the specification of variables for subprograms we can for every
argument now give its INTENT as IN, OUT or INOUT. If IN is valid, then
the actual argument can be an expression like X+Y or SIN(X) or a
constant like 37, since the value is only to be transferred to the
subprogram, but a new value is not to be returned to the calling unit.
The variables in this case may not be assigned a new value in the
subprogram. If OUT is valid, on the other hand, the actual argument has
to be a variable. At entry to the subprogram the variable is at this
stage considered to be not defined. The third case covers both
possibilities, one value on input and another on output, or possibly the
same value. Also in this case the actual argument has to be a variable.
If a variable has a pointer attribute then INTENT may not be given. The
implementation of INTENT is not yet complete in all compilers.
----------------------------------------------------------------------
- 12 -
One use for the new program unit MODULE is to take care of global
data and then it replaces the BLOCK DATA, the other use is to make a
package of new data type. As a rather large example I try to give a
package for interval arithmetic. For each value X then you have an
interval (X_lower; X_upper). At the use of the package you want to give
only the variable name X when you mean the interval. The variable X is
then supposed to be a new data type, interval. The following is on the
file interval_arithmetics.f90. Then the routine follows, written in
English.
The module continues on page 13.
MODULE INTERVAL_ARITHMETICS
TYPE INTERVAL
REAL LOWER, UPPER
END TYPE INTERVAL
INTERFACE OPERATOR (+)
MODULE PROCEDURE INTERVAL_ADDITION
END INTERFACE
INTERFACE OPERATOR (-)
MODULE PROCEDURE INTERVAL_SUBTRACTION
END INTERFACE
INTERFACE OPERATOR (*)
MODULE PROCEDURE INTERVAL_MULTIPLICATION
END INTERFACE
INTERFACE OPERATOR (/)
MODULE PROCEDURE INTERVAL_DIVISION
END INTERFACE
CONTAINS
FUNCTION INTERVAL_ADDITION(A, B)
TYPE(INTERVAL) A, B, INTERVAL_ADDITION
INTERVAL_ADDITION%LOWER = A%LOWER + B%LOWER
INTERVAL_ADDITION%UPPER = A%UPPER + B%UPPER
END FUNCTION INTERVAL_ADDITION
FUNCTION INTERVAL_SUBTRACTION(A, B)
TYPE (INTERVAL) A, B, INTERVAL_SUBTRACTION
INTERVAL_SUBTRACTION%LOWER = A%LOWER - B%UPPER
INTERVAL_SUBTRACTION%UPPER = A%UPPER - B%LOWER
END FUNCTION INTERVAL_SUBTRACTION
FUNCTION INTERVAL_MULTIPLICATION(A, B)
! POSITIVE NUMBERS ASSUMED
TYPE (INTERVAL) A, B, INTERVAL_MULTIPLICATION
INTERVAL_MULTIPLICATION%LOWER = &
A%LOWER * B%LOWER
INTERVAL_MULTIPLICATION%UPPER = &
A%UPPER * B%UPPER
END FUNCTION INTERVAL_MULTIPLICATION
FUNCTION INTERVAL_DIVISION(A, B)
! POSITIVE NUMBERS ASSUMED
TYPE(INTERVAL) A, B, INTERVAL_DIVISION
INTERVAL_DIVISION%LOWER = A%LOWER / B%UPPER
INTERVAL_DIVISION%UPPER = A%UPPER / B%LOWER
END FUNCTION INTERVAL_DIVISION
END MODULE INTERVAL_ARITHMETICS
----------------------------------------------------------------------
- 13 -
At the compilation of the above the file interval_arithmetics.mod is
created, it includes an interesting modified version of the code above.
A program that wishes to use this package includes the statement USE
INTERVAL_ARITHMETICS first among the specification statement, then the
data type INTERVAL and the four arithmetic calculations on this type are
directly available. In some cases it is desirable to only include some
of the facilities in a module. In this case you use the ONLY attribute
within the new USE statement.
USE module_name, ONLY : list_of_chosen_routines
The following is an example of a very simple main program for the
test of interval arithmetics. It is from the file interval.f90.
USE INTERVAL_ARITHMETICS
IMPLICIT NONE
TYPE (INTERVAL) : : A, B, C, D, E, F
A%LOWER = 6.9 ; A%UPPER = 7.1
B%LOWER = 10.9 ; B%UPPER = 11.1
WRITE (*,*) A, B
C = A + B ; D = A - B
E = A * B ; F = A / B
WRITE (*,*) C, D ; WRITE (*,*) E, F
END
Running this program on the Sun-computer beethoven follows:
beethoven 1 % f90 interval.f90 interval_arithmetics.f90
interval.f90 interval_arithmetics.f90:
beethoven 2 % a.out
6.9000001 7.0999999 10.8999996 11.1000004
17.7999992 18.2000008 -4.2000003 -3.7999997
75.2099991 78.8100052 0.6216216 0.6513762
beethoven 3 % exit
We compiled with the compiler f90, and the executable program was
automatically named a.out.
In a module some concepts can be defined as PRIVATE, which means that
the program units outside of this module are not able to use this
concept. Sometimes an explicit PUBLIC declaration is used, normally
PUBLIC is default. Giving the following statements
PRIVATE
PUBLIC :: VAR1
----------------------------------------------------------------------
- 14 -
it follows that all variables except VAR1 are local while VAR1 is
global. Note that both these concepts (PUBLIC and PRIVATE) either can
be given as a statement, for example
INTEGER :: IVAR
PRIVATE :: IVAR
or as an attribute
INTEGER, PRIVATE :: IVAR
and the corresponding for PUBLIC.
Exercises
---------
(7.1) Complement the module so that the package can manage
arbitrary signs on the numbers also with multiplication and
division.
(7.2) Complement the modules so that the package makes a
suitable error message when it
divides with an interval that contains zero.
(7.3) Complement also so that the local rounding error at the
operation is handled. ( This is not the case at the moment.)
8. Keyword arguments and default arguments
==========================================
Routines can now be called with keyword arguments and can use default
arguments, that means that some arguments can be given with keywords
instead of position and some do not have to be given at all, in which
case a standard or default value is used.
The use of keyword and default argument is not just as simple as it
appears from Appendix 3, section 6, program units. It is one of the
cases where an explicit interface is required. Therefore we give a
complete example here. As keywords the formal parameters of the
interface are being used, they do not have to be the same names as in
the actual subprogram. They shall not be specified in the calling
program unit.
IMPLICIT NONE
INTERFACE
SUBROUTINE SOLVE (A, B, N)
INTEGER, INTENT (IN) : : N
REAL, INTENT(OUT) : : A
REAL, INTENT(IN), OPTIONAL : : B
END SUBROUTINE SOLVE
END INTERFACE
REAL X
! Note that A, B and N are not specified as REAL or INTEGER
! in this unit.
CALL SOLVE(B=10.0,N=50,A=X)
WRITE(*,*) X
CALL SOLVE(B=10.0,N=100,A=X)
WRITE(*,*) X
CALL SOLVE(N=100,A=X)
WRITE(*,*) X
END
----------------------------------------------------------------------
- 15 -
SUBROUTINE SOLVE(A, B, N)
REAL, OPTIONAL, INTENT (IN) : : B
IF (PRESENT(B)) THEN
TEMP_B = B
ELSE
TEMP_B = 20.0
END IF
A = TEMP_B + N
RETURN
END
Note that the statement IMPLICIT NONE in the main program does not
work automatically also in the subroutine SOLVE. This subroutine
therefore ought to be given its own IMPLICIT NONE satement and the
specification of all variables used (A, B, N and TEMP_B).
It is run with the following statements
beethoven 2 % f90 program.f90
beethoven 3 % a.out
60.0000000
1.1000000E+02
1.2000000E+02
beethoven 4 %
In the last case the default argument is used, since B is not given
explicitly, it means that the default value 20 is added to the actual
argument N = 100, which gives A = 120.
It is convenient to place the interface INTERFACE in a module so the
user does not have to care so much about it. The interface is a natural
complement to the routine library. Fortran 90 looks automatically for
modules in the present directory, in the directories given in the I-list
and also in /usr/local/lib/f90 (this is the standard library for Fortran
90 using UNIX). The concept I-list is explained in Appendix 6. If you
forget INTERFACE or have an incorrect interface, usually compilation or
execution gives the error message "Segmentation error", and nothing
more.
Note that if an output variable is given as OPTIONAL and INTENT(OUT)
then you have to have it included in the argument list, if the program
at execution assigns a value to this variable. You can therefore not
use OPTIONAL in order to get a choice whether you wish to have a certain
variable outputted or not. The solution to this problem will be
discussed a little later, the PRESENT statement.
Exercises
---------
(8.1) Write a routine for the calculation of an integral of a
function. You will use keyword arguments and default
arguments so that
a) if there is no left integration limit A use the value zero.
b) if there is no right integration limit B use the value one.
c) if there is no tolerance keyword TOL use the value 0.001
for the absolute error.
(8.2) Write the interface that is required in the calling routine
in order to use the above integration routine.
----------------------------------------------------------------------
- 16 -
9. Recursion
============
A completely new possibility of Fortran 90 is recursion. Note that
it requires that you give a new property RESULT for the output variable
in the function declaration. This output variable is required inside
the functions as the "old" function name in order to store the value of
the function. At the actual call of the function, both externally and
internally, you use the outer or "old" function name. The user can
therefore ignore the output variable.
Here follows two examples, first recursive calculation of the
factorial, then recursive calculations of the Fibonacci-numbers. The
later is very inefficient. An efficient but not recursive method is
given by Brainerd, Goldberg and Adams (1990), page 226.
The listings of the routines follow, the output variables are called
FAC_RESULT and FIBO_RESULT, respectively.
RECURSIVE FUNCTION FACTORIAL(N) RESULT (FAC_RESULT)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N
INTEGER :: FAC_RESULT
IF ( N <=1 ) THEN
FAC_RESULT = 1
ELSE
FAC_RESULT = N * FACTORIAL(N-1)
END IF
END FUNCTION FACTORIAL
RECURSIVE FUNCTION FIBONACCI(N) RESULT (FIBO_RESULT)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N
INTEGER :: FIBO_RESULT
IF ( N <= 2 ) THEN
FIBO_RESULT = 1
ELSE
FIBO_RESULT = FIBONACCI(N-1) + FIBONACCI(N-2)
END IF
END FUNCTION FIBONACCI
The reason that the above calculation of the Fibonacci-numbers is so
inefficient is that each call with a certain value of N generates two
calls for the routine, which in its turn generates four calls, and so
on. Old values (or calls) are not re-used.
An interesting use of recursive technique is given by Brainerd,
Goldberg and Adams (1990), page 222, for the calculation of an
exponential function of a matrix. They give the immediate (straight
forward) expression, with the successive multiplication with a matrix,
and also a recursive variant, which can pick out the suitable squares to
optimize the calculation. Recursion is also excellent to code an
adaptive algorithm, see exercise (9.2) below.
Another very important usage of the RESULT property and the output
variable is at array valued functions. It is very easy to specify an
output variable so that it can store all the values of such a function.
Actually, it is the combination of recursive functions and array valued
functions that have forced the committee to introduce the RESULT
property.
Not only functions but also subroutines can be recursive.
----------------------------------------------------------------------
- 17 -
Exercises
---------
(9.1) Write a routine for the calculation of Tribonacci numbers.
These are formed like the Fibonacci-numbers but you start
from three numbers (all 1 at the start) and every time add
the three latest numbers to get the next one. Run
and calculate TRIBONACCI (15). Note that the calculation time
grows very fast with the argument.
(9.2) Write an adaptive routine for quadrature, i.e. calculation of
a definite integral on a certain interval.
10. Generic routines
====================
From Fortran 77 (but not from Fortran 66) we are used to have that
the elementary functions are generic, which means that a call SIN (1.0)
returns a value of type REAL, but SIN (1.0D0) returns a value with a
higher precision and of type DOUBLE PRECISION. We now also have the
possibility to write our own generic functions or subroutines. Here we
first give a complete example of a routine SWAP(A, B), which swaps the
values of variables A and B (replaces the value with each other), using
different underlying routines depending on the type of the variables
REAL, INTEGER or CHARACTER.
PROGRAM SWAP_MAIN
IMPLICIT NONE
INTEGER :: I, J, K, L
REAL :: A, B, X, Y
CHARACTER :: C, D, E, F
INTERFACE SWAP
SUBROUTINE SWAP_R(A, B)
REAL, INTENT (INOUT) :: A, B
END SUBROUTINE SWAP_R
SUBROUTINE SWAP_I(A, B)
INTEGER, INTENT (INOUT) :: A, B
END SUBROUTINE SWAP_I
SUBROUTINE SWAP_C(A, B)
CHARACTER, INTENT (INOUT) :: A, B
END SUBROUTINE SWAP_C
END INTERFACE
I = 1 ; J = 2 ; K = 100 ; L = 200
A = 7.1 ; B = 10.9 ; X = 11.1; Y = 17.0
C = 'a' ; D = 'b' ; E = '1' ; F = '"'
WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F
CALL SWAP(I, J) ; CALL SWAP(K, L)
CALL SWAP(A, B) ; CALL SWAP(X, Y)
CALL SWAP(C, D) ; CALL SWAP(E, F)
WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F
END
----------------------------------------------------------------------
- 18 -
SUBROUTINE SWAP_R(A, B)
IMPLICIT NONE
REAL, INTENT (INOUT) :: A, B
REAL :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_R
SUBROUTINE SWAP_I(A, B)
IMPLICIT NONE
INTEGER, INTENT (INOUT) :: A, B
INTEGER :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_I
SUBROUTINE SWAP_C(A, B)
IMPLICIT NONE
CHARACTER, INTENT (INOUT) :: A, B
CHARACTER :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_C
The above works very well, but the user is not so happy to have to
care with all the information about SWAP in these three different
variants in the program. The solution to this is to move everything
that has to do with the SWAP into a module and then can the module can
be used from the main program with the statement USE module name.
Please note that in the INTERFACE of the module the specific statement
MODULE PROCEDURE has to be used in order to avoid that the routines are
specified both in the INTERFACE and in the CONTAINS part. At the use
you will have to link both the module and the main program together,
e.g. with the statement
f90 part1.f90 part2.f90
Here the modules follow, it could be in the file part2.f90,
MODULE BO
INTERFACE SWAP
MODULE PROCEDURE SWAP_R, SWAP_I, SWAP_C
END INTERFACE
CONTAINS
SUBROUTINE SWAP_R(A, B)
IMPLICIT NONE
REAL, INTENT (INOUT) :: A, B
REAL :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_R
SUBROUTINE SWAP_I(A, B)
IMPLICIT NONE
INTEGER, INTENT (INOUT) :: A, B
INTEGER :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_I
----------------------------------------------------------------------
- 19 -
SUBROUTINE SWAP_C(A, B)
IMPLICIT NONE
CHARACTER, INTENT (INOUT) :: A, B
CHARACTER :: TEMP
TEMP = A ; A = B ; B = TEMP
END SUBROUTINE SWAP_C
END MODULE BO
Here follows the main program, which is now cleaned of all
uninteresting information about SWAP. It could be in the file
part1.f90.
PROGRAM SWAP_MAIN
USE BO
IMPLICIT NONE
INTEGER :: I, J, K, L
REAL :: A, B, X, Y
CHARACTER :: C, D, E, F
I = 1 ; J = 2 ; K = 100 ; L = 200
A = 7.1 ; B = 10.9 ; X = 11.1; Y = 17.0
C = ' a' ; d = 'b' ; E = '1' ; F = '"'
WRITE (*,*) I, J, K, L, A, B, C, D, E, F
CALL SWAP (I, J) ; CALL SWAP (K, L)
CALL SWAP (A, B) ; CALL SWAP (X, Y)
CALL SWAP (C, D) ; CALL SWAP (E, F)
WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F
END
11. Use of arrays and array sections
====================================
The English word "array" is translated into Swedish as "falt" or
retranslated into English as "field". We may therefore perhaps use the
word field either by mistake or as a suitable name of a specific array.
A new feature of Fortran 90 is that you can work directly with a whole
array or an array section without explicit (or implicit) DO-loops. In
the old Fortran you could in some circumstances work directly with an
whole array, but then only during I/O processing.
An array is defined to have a shape given by its number of dimensions
(called "rank") and extent for each one of these. Two arrays agree if
they have the same shape. Operations are normally done element for
element. Please note that the rank of an array is the number of
dimensions and has nothing to do with the mathematical rank of a matrix!
In the following simple example I show how you can assign matrices
with simple statements like B = A, how you can use the intrinsic matrix
multiplication MATMUL and the addition SUM and how you can use the array
sections (in the example below I use array sections who are vectors).
----------------------------------------------------------------------
- 20 -
PROGRAM ARRAY_EXAMPLE
IMPLICIT NONE
INTEGER :: I, J
REAL, DIMENSION (4,4) :: A, B, C, D, E
DO 1 = 1, 4 ! calculate a test matrix
DO J = 1, 4
A(I, J) = (I-1.2)**J
END DO
END DO
B = A*A ! element for element multiplication
CALL PRINTF(A,4) ; CALL PRINTF(B,4)
C = MATMUL(A, B) ! internal matrix multiplication
DO L = 1, 4 ! explicit matrix multiplication
DO J = 1, 4
D(I, J) = SUM( A(I,:)*B(:,J) )
END DO
END DO
CALL PRINTF(C,4) ; CALL PRINTF(D,4)
E = C - D ! comparison of the two methods
CALL PRINTF(E,4)
CONTAINS
SUBROUTINE PRINTF(A, N) ! print an array
IMPLICIT NONE
INTEGER :: N, I
REAL, DIMENSION (N, N) :: A
DO I = 1, N
WRITE(*,' (4E15.6)') A(I,:)
END DO
WRITE(*,*) ! write the blank line
END SUBROUTINE PRINTF
END PROGRAM ARRAY_EXAMPLE
As was mentioned in section 9 about recursion, functions in Fortran
90 can be array valued. In that case it is recommended to use the
RESULT property to specify the variable that is supposed to store the
array.
Fortran 90 has much larger possibilities than Fortran 77 to permit
dynamic memory allocation, which in Fortran 77 only could be done when a
sufficient storage area had been allocated in the calling program unit,
and both the array name and the required dimension(s) have to be
included as parameters in the call of the subprogram. This is the
concept adjustable array. A very simple case is the one with the last
dimension, which can be given simply with a *, assumed-size array.
Now we also have the possibilities allocatable array, automatic
array, and assumed-shape array. See the next section, Appendix 3
(section 10), and also Appendix 11 (explanation of certain terms).
Exercise
--------
(11.1) Write a routine for the solution of a system of linear
equations using Gaussian elimination with partial pivoting.
----------------------------------------------------------------------
- 21 -
12. Pointers
============
Pointers have been included in Fortran 90, but not in the usual way
as in most other languages, with pointer as a specific data type. Here
they are rather understood as an attribute to the other data types. The
reason for this new way to introduce them is that by having a special
data type for pointers, the risk of erroneous use of the pointer is very
large. A variable with a pointer attribute can be used as a usual
variable and in some new ways. Pointers in Fortran 90 are thus not
memory addresses as in other programming languages (and in certain
Fortran implementations) but rather an extra name (alias).
The increased security is obtained not only through that each
variable, which shall be used as a pointer, must be given an attribute
POINTER, but also that all variables, that will be pointed to, must be
given an attribute TARGET. An example explaining how to do this
follows.
REAL, TARGET :: B(10,10)
REAL, POINTER :: A(:,:)
A => B
The matrix B has been specified completely, i.e. with the dimensions
given explicitly. In addition, it has been stated that it can be the
target of a pointer. The matrix A, which can be used as a pointer, has
to be declared as a matrix, i.e. to be given a correct number of
dimensions, a correct rank, but the extent for this is decided later, at
the assignment (and in reality the assignment is a pointer-association)
which is done with the symbols =>. Please note that the pointer
assignment does not mean that the data in the matrix B is copied over to
the matrix A (which would have taken relatively large resources), but it
is merely a new address that is generated. To "move" data with the
pointer concept will therefore be very efficient. As an alternative the
pointer can become associated with the statement ALLOCATE, and be
deassociated with DEALLOCATE, as in the following example.
ALLOCATE (A(5,5))
DEALLOCATE (A)
There is also an internal function ASSOCIATED in order to investigate
if a pointer is associated (and if it is associated with a certain
target) and a statement NULLIFY in order to terminate the association.
IF ( ASSOCIATED (A) ) WRITE(*,*) ' A is associated '
IF ( ASSOCIATED (A, B) ) WRITE(*,*) ' A is associated with B '
NULLIFY (A)
Please remember that a pointer in Fortran 90 has both type and rank,
and that these must agree with the corresponding target. This increases
the security at the use of pointers, it is therefore not possible by
mistake to let a pointer change values of variables of other (different)
data types. The fact that you have to specify that a variable can be a
target also increases both security and efficiency of the compilation.
A simple use of pointers is to give a name to an array section.
----------------------------------------------------------------------
- 22 -
REAL, TARGET :: B(10,10)
REAL, POINTER :: A(:), C(:)
A => B(4,:) ! vector A becomes the fourth row
C => B(:,4) ! and vector C becomes the fourth
! column of the matrix B
It is not necessary to take the whole section, you can take only a
partial section. In the following example you can take a partial matrix
WINDOW of a large matrix MATRIX.
REAL, TARGET :: MATRIX(100,100)
REAL, POINTER :: WINDOW(:,:)
INTEGER :: N1, N2, M1, M2
WINDOW => MATRIX(N1:M1, N2:M2)
If you later wish to change a dimension of the partial matrix WINDOW
you only need to make a new pointer association. Please note that the
indices in WINDOW are not from N1 to M1 and from N2 to M2 but from 1 to
M1-N1+1 and from 1 to M2-N2+1.
There does not exist arrays of pointers directly in Fortran 90, but
you can construct such facilities by creating a new data type. An
example is to store a lower (or left) triangular matrix with rows with
varying length. First introduced a new data type ROW
TYPE ROW
REAL, POINTER :: R(:)
END TYPE
and then specify the two lower triangular matrices V and L as vectors
of rows with varying length
INTEGER :: N
TYPE(ROW) :: V(N), L(N)
after which you can allocate the matrix V as below (and in the
corresponding way you can allocate the matrix L)
DO I = 1, N
ALLOCATE (V(I)%R(1:I)) ! Various length of rows
END DO
The statement V = L then becomes equivalent with V(I)%R => L(I)%R for
all the components, i.e. all values of I. Please note that in this
application there is no TARGET required.
You have to be careful when you use pointers. In the following
simple example we look at ordinary scalar floating-point numbers.
----------------------------------------------------------------------
- 23 -
REAL, TARGET :: A
REAL, POINTER :: P, Q
A = 3.1416
P => A
Q => P
A = 2.718
WRITE(*,*) Q
Here the value of Q equals 2.718 since both P and Q point towards the
same variable A and that one has just changed its value from 3.1416 to
2.718. We now make a simple variation.
REAL, TARGET :: A, B
REAL, POINTER :: P, Q
A = 3.1416
B = 2.718
P => A
Q => B
Now both the values of A and P are equal to 3.1416 and the values of
both B and Q are 2.718. If we now give the statement Q = P, all four
variables will get the value 3.1416, which means that an ordinary
assignment of pointer variables has the same effect as the conventional
assignment B = A. If we instead give a pointer association Q => P, then
the three variables A, P and Q all have the value 3.1416, while B
contains the value 2.718. In the second case Q only points to the same
variable as P while in the first case Q becomes the same as P, and the
value addressed by Q becomes equal to the value addressed by P.
Important application of pointers are lists and trees, and
especially dynamic arrays.
Exercises
---------
(12.1) Use pointers in order to assign all even elements of a vector
the value 13 and all odd elements of a vector the value 17.
(12.2) Specify two pointers, and let one of them point to a whole
vector and the other one point to the seventh element of the
same vector.
(12.3) Use pointers to specify a vector in such a way, that it is
given its size (its extent) in a subroutine but can be used in
the main program. This is one implementation of dynamic memory
allocation.
----------------------------------------------------------------------
- 24 -
13. The new precision concept
==============================
The problem with the older versions of Fortran was that simple
precision on one computer could correspond to a higher precision or
DOUBLE precision on another computer and that the data type DOUBLE
PRECISION COMPLEX or COMPLEX*16 was not available in all systems (and of
course not in the standard).
In Fortran 90 there are standard functions to check the precision of
variables (see Appendix 5, section 8, where for example PRECISION(X)
gives the number of significant digits in numbers of the same kind as
the variable X). In Fortran 90 there are also possibilities to specify
for each variable how many significant digits can be used with the
floating-point numbers of this kind. The two common precisions single
precision (SP) and double precision (DP) on a system based on IEEE 754
can specified with
INTEGER, PARAMETER :: SP = SELECTED_ REAL_ KIND(6,37)
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,307)
REAL(KIND=SP) :: single_precision_variables
REAL(KIND=DP) :: double_precision_variables
If we wish to work with at least 14 decimal digits accuracy and at
least decimal exponents between - 300 and + 300 we can choose the
following integer parameters
INTEGER, PARAMETER :: WP = SELECTED_REAL_KIND(14,300)
REAL(KIND=WP) :: working_precision_variables
Regrettably now we have to give all floating point constants with the
additional suffix _WP, for example
REAL(KIND=WP) :: PI
PI = 3.141592653589793_WP
while since the intrinsic functions are generic, they will
automatically use the correct data type and kind, which means that the
argument determines which kind the result should have (usually the same
as the argument).
With this method you will in practice obtain double precision on
systems based on IEEE 754 and single precision on computers like Cray or
computers based based on the Digital Equipment Alpha-processor, which in
all cases means a precision of about 15 significant digits.
14. Additional problems at the transition
=========================================
Removal of automatic generation of new lines at input
-----------------------------------------------------
A problem, that can arise when moving from Fortran 77 to Fortran 90,
is dependent on a usual deviation from the standard, and has to do
especially with user programs. What we consider here, is the use in the
FORMAT of the dollar symbol $ in order to remove the generation of a new
line (Line Feed/Carriage Return), before the user gives a new value to a
variable. This is a common extension of Fortran 77, but it generates a
compilation error in Fortran 90 for the dollar symbol and therefore
another solution has to be found.
----------------------------------------------------------------------
- 25 -
With many Fortran 77 implementations we wrote
PROGRAM TEST
REAL X
WRITE(*,10)
10 FORMAT('Give X = ',$)
READ(*,*) X
WRITE(*,*) X
END
In Fortran 90 we use "non-advancing I/O" instead. We therefore write
the following
PROGRAM TEST
IMPLICIT NONE
REAL X
WRITE(*,'(A)',ADVANCE='NO') 'Give X = '
READ(*,*) X
WRITE(*,*) X
END
Both those programs give the same result, you may give the value of
the variable on the same row as the text "Give X = ". Non-advancing I/O
can not be used with list-directed I/O or on NAMELIST.
Varying system for the treatment of matrices
---------------------------------------------
In Fortran 77 there is no dynamic memory allocation, you therefore
have to give a sufficiently large dimension in the calling program and
keep in mind the "leading dimension" in the called program unit. Now
when you use Fortran 90 you prefer to have an array of the same shape
and size as the logical size of a matrix. An assignment that performs
this transformation is easy to achieve. We assume that a quadratic
matrix from the calling program unit is available in the subprogram
under consideration as the array A with the dimensions A(IA, *) and that
we wish to move this to an array B with the dimension specification B(N,
N), where N at the same time is the mathematical dimension of the
matrix. The following assignment statement gives what you want,
provided IA is not less than N.
B = A(1:N, 1:N)
Differences in the use of logical variables
-------------------------------------------
The new standard contains many more words and is more explicit, while
some of the rules are formulated in such a way that the probability that
they are obeyed by the compiler or rather the compiler writer is much
greater. An example is comparison of a logical variables, which under
Fortran 90 has to be done with .EQV. or .NEQV. while this in practice
often was possible in Fortran 77 also with .EQ. and .NE. If you try to
perform such a comparison in Fortran 90 with the new equality symbol ==
you can get a confusing error message, complaining that .EQ. may not be
used in this context. The reason is of course that == is just an
alternative spelling of .EQ.
----------------------------------------------------------------------
- 26 -
Small things of importance
--------------------------
It has been rather common to use the variable name SUM in order to
store the temporary value at the summation in a DO-loop. This name is
now less convenient to use, since SUM is also the name of an automatic
summation, see Appendix 5 (section 14 on array functions). Other
dangerous variable names can be ALL, HUGE, INDEX, INT, KIND, MASK,
SCALE, SIZE, TINY and TRIM. If you use a name that is being used by the
Fortran language, the normal effect is that the intrinsic function is no
longer available.
In some old Fortran dialects the statement TYPE was used to write on
a typewriter terminal and PRINT was used to write on the line printer.
The concept TYPE now has a completely new meaning in Fortran 90, to
declare user-defined data types.
15. Use of program libraries
============================
A problem with Fortran 90 is that so far not so many programs or
program libraries are available in the language. We can therefore be
forced to use "old" program libraries usually from Fortran 77. This can
be done in two ways.
Method 1. A simple method is to recompile the Fortran 77 routines
with the Fortran 90 compiler using fixed form. This should work since
Fortran 77 is a pure subset of Fortran 90. The problem is that the
library does not always exactly follow the standard. NAG recognizes
specifically that although the data type DOUBLE PRECISION COMPLEX or
COMPLEX*16 is not in the standard, it is in many implementations of
Fortran 77. A subprogram which uses this extension therefore has to be
modified. Also possible assignments of binary, octal or hexadecimal
constants is difficult since they are standardized only in Fortran 90.
Also note that the very common notation REAL*8 for DOUBLE PRECISION is
not permitted in Fortran 90.
Method 2. The second method is to link the Fortran 77 library
directly with the Fortran 90 programs. The practical problem was
previously that the NAG's Fortran 90 system f90 under version 1.1 did
not include support for linkage of libraries, therefore you first had to
compile with the Fortran compiler using f90 -c and then link with the C
compiler using cc. NAG recommended the following commands
f90 -c test.f
cc test.o -lnag -o a.out /usr/local/lib/f90/f90rt0.o\
/usr/local/lib/f90/libf90.a -lI77 -lF77 - lm
which worked on Sun after that we had removed -lI77 since we did not
find this library on the system here in Linkoping.
From version 1.2 it is much simpler because linkage support is now
included and on the Sun you only give the following command
f90 test.f -lnag - lF77
while on DEC you give the command
f90 test.f lnag -lfor -lutil -li -lots
----------------------------------------------------------------------
- 27 -
Thus a few more libraries have to be included on the DEC.
In addition to the problems mentioned above we can also note that
Fortran allows routine names as arguments, which may be treated
differently in the NAG's Fortran 90 compiler based on translation to C
and in the existing Fortran 77 compiler. There is therefore a danger
for linkage errors. On the Sun using Sun OS 4.2.1 we also get
completely wrong results when using library functions in single
precision, since the C system on Sun converts all function calls into
double precision. NAG claims to have solved this problem with single
precision in release 2.1.
The method can thus be used when neither complex variables in double
precision nor routine names are used as arguments, and provided that you
are not using single precision library functions on a Sun. It is very
essential to get all the required libraries at the linkage process,
those are for example the mathematical libraries in Fortran 77, Fortran
90 and C. Sometimes the libraries have to be given in the correct
order.
16. Peculiarities in the language Fortran 90
============================================
o Using the free source form it is not permitted to have comments
starting with a C or a * in column 1. This would be a violation
of the free format. The new character ! (exclamation mark) has
to be used.
o The repetition factor at the assignment of initial values was
permitted in the old DATA-statement, but it is not permitted in
the new direct assignment within (/ /). It is also not permitted
in connection with the PARAMETER statement or attribute.
o An advantage compared with earlier Fortran standards is that
the standard now requires that the compiler can signal if the
user deviates from the permitted standard.
It is required that a Fortran 90 compiler can signal
- use of syntax not defined in the standard
- violation of the syntax rules
- use of kinds not available
- use of obsolete constructs (or statements)
- use of non-Fortran characters (for example, Swedish or Cyrillic
characters) outside of character strings or comments
- violation of the area of validity for variable names, names of
the DO-loops and the corresponding names like IF, CASE and operators.
- the reason that a program is not accepted by the compiler
The above means that it is permitted to include extensions to
Fortran 90. It has to be possible to ask the program to signal for
any extensions outside the Fortran 90 standard.
o The compiler often performs some data flow analysis.
Answers and comments to the user exercises
==========================================
It is assumed that when nothing else is stated, the implicit rules
about integers and floating-point numbers are used, i.e. IMPLICIT NONE
has not been used. In those cases where runs on computers shall be
done, I refer the reader to Appendix 6 for MS-DOS and UNIX computers.
----------------------------------------------------------------------
- 28 -
This Appendix is however mainly about the use of the Fortran 90 system
from NAG. For IBM PC there is a very complete documentation in the
booklet "NAGware FTN90 compiler". In all other cases please try the
manuals from the computer or compiler manufacturer.
(1.1) If the compilation or the execution run fails it is probably
an error already in the Fortran 77 program, or you have used
some extension to standard Fortran 77.
(1.2) In this fails it probably depends on that some incorrect
commands were interpreted as variables when using fixed form,
but now when blanks are significant these syntax errors are
discovered.
(1.3) Fortran 77 does not give any error either on the compilation
or execution. Compilation in Fortran 90 fixed form gives a
warning from the compiler that the variable ZENDIF is used
without being assigned any value. The program is interpreted
in such a way that THENZ, ELSEY, and ZENDIF becomes
ordinary floating-point variables. Compilation in Fortran 90
free form, however, gives a number of syntax errors. The
correct version of the program shall contain three extra
carriages return as below.
LOGICAL L
L = .FALSE.
IF (L) THEN
Z = 1.0
ELSE
Y=Z
END IF
END
REMARK: Also certain Fortran 77 compilers give a warning about the
variable ZENDIF, which has not been assigned any value.
(2.1) Using fixed form it means LOGICAL L, i.e. the variable L is
specified as logical. Using free form you will get a syntax
error.
(2.2) REAL, PARAMETER :: K = 0.75
(2.3) INTEGER, DIMENSION(3,4) :: PELLE
(2.4) INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,99)
(2.5) REAL (KIND=DP) :: E, PI
(2.6) E = 2.718281828459045_DP
PI = 3.141592653589793_DP
(2.7) No, it is not correct since a comma is missing between
REAL and DIMENSION. In the form it has been written, the
statement is interpreted as a specification of the old type
of the floating-point matrix DIMENSION (with the specified
dimensions), and an implicit specification of the new type of
a scalar floating-point number AA. Formally, it is a correct
specification. The variable name DIMENSION is permitted in
Fortran 90, just as the variable name REAL is permitted in
both Fortran 77 and Fortran 90, but both should be avoided.
The variable name DIMENSION is of course too long in
standard Fortran 77.
----------------------------------------------------------------------
- 29 -
(2.8) Yes, it is correct, but it is not suitable since it kills the
intrinsic function REAL for explicit conversion of a variable
of another type to the type REAL. It is however nothing that
prevents you from using a variable of the type REAL with the
name REAL, since Fortran does not have reserved words.
(2.9) No, it is not correct, at COMMON you do not use the double
colon at the specification. The correct specification is the
old familiar one.
COMMON A
(3.1) Variables A and B are assigned the specified values, but the
whole rest of the line becomes a comment.
(3.2) No, on the second row the blank space after the ampersand (&)
is not permitted. It interrupts the identifier ATAN into two
identifiers AT and AN. If the blank is removed the two lines
become correct. Free form is assumed, since & is not a
continuation character in fixed form.
(4.1) The statement is not permitted, which however is shown not at
compilation but at execution. You can instead write
WRITE(*,*) ' HI '
or
WRITE(*,'(A)') ' HI '
which both write out the text HI on the standard unit for
output. You are not permitted to give the text, which you
wish to to print, directly where the output format is to be
given.
(4.2) They write large and small numbers with an integer digit, six
decimals and an exponent, while numbers in between are
written in the natural way. In this case we thus get
1.000000E-03
1.00000
1.000000E+06
Numbers from 0.1 to 100 000 are written in the natural way
and with six significant digits.
(6.1) SELECT CASE (N)
CASE(:-1)
! Case 1
CASE(0)
! Case 2
CASE(3,5,7,11,13)
! Case 3
END SELECT
(6.2) SUMMA = 0.0
DO I = 1, 100
IF ( X(I) == 0.0) EXIT
IF ( X(I) < 0.0) CYCLE
SUMMA = SUMMA + SQRT (X(I))
END
----------------------------------------------------------------------------
- 30 -
The English word sum is not suited as the variable name in this case,
since this is also an intrinsic function. Summa is the Swedish word for
sum.
(7.1) Use the functions MIN and MAX to find the largest
and smallest values of all the combinations
A%LOWER * B%LOWER, A%LOWER * B%UPPER, A%UPPER * B%LOWER and
A%UPPER * B%UPPER
at multiplication and the corresponding at division.
(7.2) Test if B%LOWER <= 0 <= B%UPPER in which case an error
message shall be given.
(7.3) Since we do not have direct access to machine arithmetics
(i.e. commands of the type round down or round up) you can
get a reasonable simulation through subtraction and addition
with the rounding constant. In principle the effect of
rounding is then doubled.
(8.1) SUBROUTINE SOLVE(F, A, B, TOL, EST, RESULT)
REAL, EXTERNAL :: F
REAL, OPTIONAL, INTENT (IN) :: A
REAL, OPTIONAL, INTENT (IN) :: B
REAL, OPTIONAL, INTENT (IN) :: TOL
REAL, INTENT(OUT), OPTIONAL :: EST
REAL, INTENT(OUT) :: RESULT
IF (PRESENT(A)) THEN
TEMP_A = A
ELSE
TEMP_A = 0.0
END IF
IF (PRESENT(B)) THEN
TEMP_B = B
ELSE
TEMP_B = 1.0
END IF
IF (PRESENT(TOL)) THEN
TEMP_TOL = TOL
ELSE
TEMP_TOL = 0.001
END IF
! Here the real calculation should be, but it is here replaced
! with the simplest possible approximation, namely the middle
! point approximation without an error estimate.
RESULT = (TEMP_B - TEMP_A)&
* F(0.5*(TEMP_A+TEMP_B))
IF (PRESENT(EST)) EST = TEMP_TOL
RETURN
END SUBROUTINE SOLVE
The very simple integral calculation above can be replaced
by the adaptive quadrature in exercise (9.2).
----------------------------------------------------------------------
- 31 -
(8.2) INTERFACE
SUBROUTINE SOLVE (F, A, B, TOL, EST, RESULT)
REAL, EXTERNAL :: F
REAL, INTENT(IN), OPTIONAL :: A
REAL, INTENT(IN), OPTIONAL :: B
REAL, INTENT(IN), OPTIONAL :: TOL
REAL, INTENT(OUT), OPTIONAL :: EST
REAL, INTENT(OUT) :: RESULT
END SUBROUTINE SOLVE
END INTERFACE
(9.1) RECURSIVE FUNCTION TRIBONACCI (N) RESULT (T_RESULT)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N
INTEGER :: T_RESULT
IF ( N <= 3 ) THEN
T_RESULT = 1
ELSE
T_RESULT = TRIBONACCI(N-1 )+ &
TRIBONACCI(N-2) + TRIBONACCI(N-3)
END IF
END FUNCTION TRIBONACCI
The calling program or main program can be written
IMPLICIT NONE
INTEGER :: N, M, TRIBONACCI
N = 1
DO
IF ( N <= 0 ) EXIT
WRITE (*,*) ' GIVE N '
READ(*,*) N
M = TRIBONACCI (N)
WRITE(*,*) N, M
END DO
END
and gives the result TRIBONACCI(15) = 2209.
(9.2) The file quad.f90 below contains a function for adaptive
numerical quadrature (integration). We use the trapezoidal formula,
divide the step size with two, and perform Richardson extrapolation.
The method is therefore equivalent to the Simpson formula. As an error
estimate we use the model in Linkoping, where the error is assumed less
than the modulus of the difference between the two not extrapolated
values. If the estimated error is too large, the routine is applied
once again on each of the two subintervals, in that case the permitted
error in each one of the subintervals becomes half of the error
previously used.
----------------------------------------------------------------------
- 32 -
RECURSIVE FUNCTION ADAPTIVE_QUAD (F, A, B, TOL, ABS_ERROR) &
RESULT (RESULT)
IMPLICIT NONE
INTERFACE
FUNCTION F(X) RESULT (FUNCTION_VALUE)
REAL, INTENT(IN) :: X
REAL :: FUNCTION_VALUE
END FUNCTION F
END INTERFACE
REAL, INTENT(IN) :: A, B, TOL
REAL, INTENT(OUT) :: ABS_ERROR
REAL :: RESULT
REAL :: STEP, MIDDLE_POINT
REAL :: ONE_TRAPEZOIDAL_AREA, TWO_TRAPEZOIDAL_AREAS
REAL :: LEFT_AREA, RIGHT_AREA
REAL :: DIFF, ABS_ERROR_L, ABS_ERROR_R
STEP = B-A
MIDDLE_POINT= 0.5 * (A+B)
ONE_TRAPEZOIDAL_AREA = STEP * 0.5 * (F(A)+ F(B))
TWO_TRAPEZOIDAL_AREAS = STEP * 0.25 * (F(A) + F(MIDDLE_POINT))+&
STEP * 0.25 * (F(MIDDLE_POINT) + F(B))
DIFF = TWO_TRAPEZOIDAL_AREAS - ONE_TRAPEZOIDAL_AREA
IF ( ABS (DIFF) < TOL ) THEN
RESULT = TWO_TRAPEZOIDAL_AREAS + DIFF/3.0
ABS_ERROR = ABS(DIFF)
ELSE
LEFT_AREA = ADAPTIVE_QUAD (F, A, MIDDLE_POINT, &
0.5*TOL, ABS_ERROR_L)
RIGHT_AREA = ADAPTIVE_QUAD (F, MIDDLE_POINT, B, &
0.5*TOL, ABS_ERROR_R)
RESULT = LEFT_AREA + RIGHT_AREA
ABS_ERROR = ABS_ERROR_L + ABS_ERROR_R
END IF
END FUNCTION ADAPTIVE_QUAD
The file test_quad.f90 for the test of the above routine for adaptive
numerical quadrature requires an INTERFACE both for the function F and
for the quadrature routine ADAPTIVE_QUAD. Note that for the latter you
must specify the function both REAL and EXTERNAL and that routine
follows on the next page, page 33.
----------------------------------------------------------------------
- 33 -
PROGRAM TEST_ADAPTIVE_QUAD
IMPLICIT NONE
INTERFACE
FUNCTION F(X) RESULT (FUNCTION_VALUE)
REAL, INTENT(IN) :: X
REAL :: FUNCTION_VALUE
END FUNCTION F
END INTERFACE
INTERFACE
RECURSIVE FUNCTION ADAPTIVE_QUAD &
(F, A, B, TOL, ABS_ERROR) RESULT (RESULT)
REAL, EXTERNAL :: F
REAL, INTENT (IN) :: A, B, TOL
REAL, INTENT (OUT) :: ABS_ERROR
REAL :: RESULT
END FUNCTION ADAPTIVE_QUAD
END INTERFACE
REAL :: A, B, TOL
REAL :: ABS_ERROR
REAL :: RESULT, PI
INTEGER :: I
PI = 4.0 * ATAN(1.0)
A= -5.0
B = +5.0
TOL =0.1
DO I = 1, 5
TOL = TOL/10.0
RESULT = ADAPTIVE_QUAD (F, A, B, TOL, ABS_ERROR)
WRITE(*,*)
WRITE(*,"(A, F15.10, A, F15.10)") &
"The integral is approximately ", &
RESULT, "with approximate error estimate ", &
ABS_ERROR
WRITE(*,"(A, F15.10, A, F15.10)") &
"The integral is more exactly ", &
SQRT(PI), " with real error ", &
RESULT - SQRT(PI)
END DO
END PROGRAM TEST_ADAPTIVE_QUAD
We are of course not permitted to forget the integrand, which we
prefer to put in the same file as the main program. Declarations are of
the new type especially with respect to that the result is returned in a
special variable.
----------------------------------------------------------------------
- 34 -
FUNCTION F(X) RESULT (FUNCTION_VALUE)
IMPLICIT NONE
REAL, INTENT(IN) :: X
REAL :: FUNCTION_VALUE
FUNCTION_VALUE = EXP(-X**2)
END FUNCTION F
Now it is time to do the test on the Sun computer named beethoven. I
have adapted the output a little in order to get it more compact. The
error estimated is rather realistic, at least with this integrand, which
is the unnormalized error function.
beethoven 3 % f90 test_quad.f90 quad.f90
test_quad.f90:
quad.f90:
beethoven 4 % a.out
The integral is 1.7733453512 with error estimate 0.0049186843
The integral is 1.7724539042 with real error 0.0008914471
The integral is 1.7724548578 with error estimate 0.0003375171
The integral is 1.7724539042 with real error 0.0000009537
The integral is 1.7724541426 with error estimate 0.0000356939
The integral is 1.7724539042 with real error 0.0000002384
The integral is 1.7724540234 with error estimate 0.0000046571
The integral is 1.7724539042 with real error 0.0000001192
The integral is 1.7724539042 with error estimate 0.0000004876
The integral is 1.7724539042 with real error 0.0000000000
beethoven 5 %
In the specification above of the RECURSIVE FUNCTION ADAPTIVE_QUAD
you may replace the line
REAL, EXTERNAL :: F
with a complete repetition of the interface for the integrand function,
INTERFACE
FUNCTION F(X) RESULT (FUNCTION_VALUE)
REAL, INTENT(IN) :: X
REAL :: FUNCTION_VALUE
END FUNCTION F
END INTERFACE
With this method an explicit EXTERNAL statement is no longer
required, but you get a nested INTERFACE.
(11.1)
SUBROUTINE SOLVE_SYSTEM_OF_LINEAR_EQUATIONS(A, X, B, ERROR)
IMPLICIT NONE
! Array specifications
REAL, DIMENSION (:, :), INTENT (IN) :: A
REAL, DIMENSION (:), INTENT (OUT):: X
REAL, DIMENSION (:), INTENT (IN) :: B
LOGICAL, INTENT (OUT) :: ERROR
! The working area M is A expanded with B
REAL, DIMENSION (SIZE (B), SIZE (B) + 1) :: M
INTEGER, DIMENSION (1) :: MAX_LOC
REAL, DIMENSION (SIZE (B) + 1) :: TEMP_ROW
INTEGER :: N, K, I
----------------------------------------------------------------------
- 35 -
! Initializing M
N = SIZE (B)
M (1:N, 1:N) = A
M (1:N, N+1) = B
! Triangularization
ERROR = .FALSE.
TRIANGULARIZATION_LOOP: DO K = 1, N - 1
! Pivoting
MAX_LOC = MAXLOC (ABS (M (K:N, K)))
IF ( MAX_LOC(1) /= 1 ) THEN
TEMP_ROW (K:N+1 ) =M (K, K:N+1)
M (K, K:N+1)= M (K-1+MAX_LOC(1), K:N+1)
M (K-1+MAX_LOG(1), K:N+1) = TEMP_ROW(K:N+1)
END IF
IF (M (K, K) == 0) THEN
ERROR = .TRUE. ! Singular matrix A
EXIT TRIANGULARIZATION_LOOP
ELSE
TEMP_ROW (K+1:N) = M (K+1:N, K) / M (K, K)
DO I = K+1, N
M (I, K+1:N+1) = M (I, K+1:N+1) - &
TEMP_ROW (I) * M (K, K+1:N+1)
END DO
M (K+1:N, K) =0 ! These values are not used
END IF
END DO TRIANGULARIZATION_LOOP
IF (M, (N, N) == 0) ERROR = .TRUE. ! Singular matrix A
! Re-substitution
IF (ERROR) THEN
X = 0.0
ELSE
DO K = N, 1, -1
X (K) = (M (K, N+1) - &
SUM (M (K, K+1:N)* X (K+1:N)) ) / M (K, K)
END DO
END IF
END SUBROUTINE SOLVE_SYSTEM_OF_LINEAR_EQUATIONS
The input matrix A and the vectors B and X are specified as
assumed-shape arrays, i.e. type, rank and name are given here, while
the extent is given in the calling program unit, using an explicit
interface.
Please note that the intrinsic function MAXLOC as a result gives an
integer vector, with the same number of elements as the rank (number of
dimensions) of the argument. In our case the argument is a vector and
therefore the rank is 1 and MAXLOC is a vector with only 1 element.
This is the reason why the local variable MAX_LOC has been declared as a
vector with 1 element. If you declare MAX_LOC as a scalar you get a
compilation error. The value of course is the index for the largest
element (the first of the largest if there are several of these).
----------------------------------------------------------------------
- 36 -
Also note that the numbering starts with 1, in spite of that we are
looking at the vector with the elements running from K to N. I prefer
not to perform the pivoting process (that is the actual exchange of
rows) in the special case that the routine finds that the rows already
are correctly located, i.e. when MAX_LOC(1) is 1.
The calculation is interrupted as soon as a singularity is found.
Please note that this can occur so late that it is not noted inside the
loop, thus the extra check immediately after the loop, for the final
element M(N, N).
At the pivoting process we use the vector TEMP_ROW first at the
exchange of lines, then also to store the multipliers in the Gauss
elimination.
In this first version we only use array sections of vector type at
the calculations, but we will now introduce the function SPREAD in order
to use array sections of matrix type, and in this case the explicit
inner loop disappears (DO I = K+1, N).
The function SPREAD(SOURCE, DIM, NCOPIES) returns an array of the
same type as the argument SOURCE, but with the rank increased by one.
The parameters DIM and NCOPIES are integers. If NCOPIES is negative the
value zero is used instead. If SOURCE is a scalar, SPREAD becomes a
vector with NCOPIES elements which all have the same value as SOURCE.
The parameter DIM gives which index is to be increased. That must be a
value between 1 and 1+(rank of SOURCE). If SOURCE is a scalar, then DIM
has to be 1. The parameter NCOPIES gives the total number of elements
in the new dimension, thus not only the number of new copies but also
the original.
If the variable A corresponds to the following array (/ 2, 3, 4 /) we get
SPREAD (A, DIM=1, NCOPIES=3) SPREAD (A, DIM=2, NCOPIES=3)
2 3 4 2 2 2
2 3 4 3 3 3
2 3 4 4 4 4
I now use array sections of matrix type through replacing the inner loop,
DO I = K+1, N
M (I, K+1:N+1) = M (I, K+1:N+1) - &
TEMP_ROW (I) * M (K, K+1:N+1)
END DO
with
M (K+1:N, K+1:N+1) = M (K+1:N, K+1:N+1) &
- SPREAD( TEMP_ROW (K+1:N), 2, N-K+1) &
* SPREAD( M (K, K+1:N+1), 1, N-K)
The reason that we have to make it almost into a muddle with the
function SPREAD is that in the explicit loop (at a fixed value of I) the
variable TEMP_ROW(I) is a scalar constant, which is multiplied by N-K+1
different elements of the matrix M, or a vector of M. On the other
hand, the same vector of M is used for all N-K values of I. The
rearrangement of the matrices has to be done to obtain two matrices with
the same shape as the submatrix M( K+1:N, K+1:N+1), that is N-K rows and
N-K+1 columns, since all calculations on arrays in Fortran 90 are
element by element.
----------------------------------------------------------------------
- 37 -
Unfortunately it is rather difficult to get the parameters to the
intrinsic function SPREAD absolutely correct. In order to get them
correct you can utilize the functions LBOUND and UBOUND in order to
obtain the lower and upper dimension limits, and you can also write the
new array with the following statement
WRITE(*,*) SPREAD (SOURCE, DIM, NCOPIES)
Please note that the output is done column by column (i.e. the first
index is varying fastest, as it is usual in Fortran). You can use the
lower and upper dimension limits for more explicit output statements
that give an output which is better suited to how the array looks.
However, here you have to first make an assignment to an array,
specified in the usual way with the correct shape, in order to use the
indices in the ordinary way. Please remember that the indices in a
construct like SPREAD automatically go from one as the lower limit.
Even when you give something like A(4:7) as SOURCE then the result will
have the index going or ranging from 1 to 4.
(12.1) We assume that the vector has a fixed dimension, and we
perform a control output of a few of the values.
REAL, TARGET, DIMENSION(1:100) :: VECTOR
REAL, POINTER, DIMENSION(:) :: ODD, EVEN
ODD => VECTOR(1:100:2)
EVEN => VECTOR(2:100:2)
EVEN = 13
ODD = 17
WRITE(*,*) VECTOR(11), VECTOR(64)
END
(12.2) We assume that the given vector has a fixed dimension.
REAL, TARGET, DIMENSION(1:10) :: VECTOR
REAL, POINTER, DIMENSION(:) :: POINTER1
REAL, POINTER :: POINTER2
POINTER1 => VECTOR
POINTER2 => VECTOR(7)
(12.3) We use an INTERFACE with pointers in the main program and
allocate, using pointers, a vector in the subroutine. In
this way we get a dynamically allocated vector.
PROGRAM MAIN_PROGRAM
INTERFACE
SUBROUTINE SUB(B)
REAL, DIMENSION (:), POINTER :: B
END SUBROUTINE SUB
END INTERFACE
REAL, DIMENSION (:), POINTER :: A
CALL SUB(A)
! Now we can use the vector A.
! Its dimension was determined in the subroutine,
! the number of elements is available as SIZE(A).
END PROGRAM MAIN_PROGRAM
----------------------------------------------------------------------
- 38 -
SUBROUTINE SUB(B)
REAL, DIMENSION (:), POINTER :: B
INTEGER M
! Now we can assign a value to M, for example through
! an input statement.
! When M has been assigned we can allocate B as a vector.
ALLOCATE (B(M))
! Now we can use the vector B.
END SUBROUTINE SUB
Note: The method above is even more useful for allocating
matrices.
References
==========
The comments below are by Bo Einarsson.
1. Jeanne C. Adams, Walter S. Brainerd, Jeanne T. Martin, Brian T.
Smith and Jerrold L. Wagener: Fortran 90 Handbook, Complete
ANSI/ISO Reference, McGraw-Hill, New York 1992. ISBN 0-07-000406-4.
$79.95.
Complete guide to Fortran 90 and its use. Written by persons that
were involved in the development of Fortran 90. Contains hundreds
of examples. However, most of these are very short and not complete
program units. Much more readable and easier to use than
the formal standards, but in spite of this it is not suitable as the
only aid to a beginner in Fortran 90.
2. ANSI: Programming Language Fortran, X3.9-1978, American National
Standard. $24.00.
The official standard for Fortran 77. It is possible to use as for
reference, but it requires that you know the basics of the language.
3. ANSI: Programming Language Fortran 90, X3.198-1992, American
National Standard.
The official American standard for Fortran 90. The same book as
ISO below (reference 14).
4. Walter S. Brainerd, Charles H. Goldberg and Jeanne C. Adams:
Programmer's Guide to Fortran 90, Second Edition, Unicomp, Albuquerque,
New Mexico, 1994. $30.00.
One of the first books about Fortran 90. Easy to read. Each
new concept that is presented is given a simple example
and therefore you can easily see how each concept is used.
It also treats the old Fortran 77 in the last chapter. The book
is written by persons that were involved in the development
of Fortran 90. The book is recommended.
5. Thomas F. Coleman and Charles Van Loan: Handbook for Matrix
Computations, Frontiers in Applied Mathematics, Vol. 4, SIAM,
Philadelphia 1988. ISBN 0-89871-227-0.
The first chapter is an excellent introduction to Fortran 77. Very
easy to read. Also treats BLAS, LINPACK and MATLAB.
6. Martin Counihan: Fortran 90, Pitman, London 1990. ISBN 0-273-03073-6.
I have not seen this book, but it is rumoured to be easy to understand
and it gives a lot of examples.
----------------------------------------------------------------------
- 39 -
7. Cray: CF77 Compiling System, Volume 1: Fortran Reference Manual,
SR-3071 5.0, June 1991.
Very thick book. Treats not only the whole language Fortran 77 but
also how it is used on the Cray. In the same series are manuals
about error messages, vectorization and parallelization.
8. Cray: CF90 Fortran Language Reference Manual, SR-3902 1.0, December 1993.
Very thick book. Treats not only the whole language Fortran 90 but
also how it is used on the Cray, with some extensions.
9. DEC: DEC Fortran, Language Reference Manual, AA-PNU0A-TK, March
1992.
This is a complete manual which also treats Fortran 77 and all
the extensions made by Digital. Necessary for DEC-programmers.
Very expensive.
10. DEC: DEC Fortran for ULTRIX RISC Systems, User Manual,
AA-PNU1A-TE, March 1992.
Auxiliary manual on the ULTRIX - environment for Fortran 77. Necessary
book for the serious DEC-programmer. Is usually bought together
with the book above.
11. Stacey L. Edgar: FORTRAN for the '90s, Problem Solving for
Scientists and Engineers, Computer Science Press, New York, 1992.
ISBN 0-7167-8247-2. $19.95.
Complete textbook in both programming in Fortran 77 and in Fortran
90. Many examples from many different areas from science and
technology. In each chapter new features of Fortran 90 are discussed
and Fortran 90 is also more fully discussed in the concluding
chapter. The book is recommended.
12. Torgil Ekman and Goran Eriksson: Programmering i Fortran 77,
Third edition, Studentlitteratur, Lund 1984. ISBN 91-44-16663-X
An excellent tutorial on Fortran 77. Describes all the
commands of Fortran. It is recommended to previously have read
a book on another language, like Pascal. Appendix C is both
well-done and very important. The book is recommended to those
who are fluent in Swedish.
13. T. M. R. Ellis: Fortran 77 Programming, Second Edition,
Addison-Wesley Publishing Company, Reading, Massachusetts 1990.
ISBN 0-201-41638-7.
Complete book in both programming in general and in Fortran 77.
Many examples and good exercises. The last chapter treats
Fortran 90.
14. T. M. R. Ellis, I. R. Philips and T. M. Lahey: Fortran 90
Programming, Addison-Wesley Publishing Company, Reading,
Massachusetts 1994. ISBN 0-201-54446-6.
Complete book in both programming in general and in Fortran 90.
Many examples and good exercises. The book is recommended.
15. High Performance Fortran Forum: High Performance Fortran
Language Specification, Version 1.0, 3 May 1993. Technical Report
CRPC-TR 92225, Center for Research on Parallel Computation, Rice
University, Houston, Texas 77251.
Available via anonymous ftp from titan.cs.rice.edu as the file
/public/HPFF/draft/hpf-v10-final.ps.Z. Includes 12 + 184 pages.
Very easy to read compared with most other standards, and has many
good examples.
----------------------------------------------------------------------
- 40 -
16. Wilhelm Gehrke: Fortran 90 Referenz-Handbuch, Hanser, Munich
1991. ISBN 3-446-16321-2. DM 168.00.
Complete description in German of Fortran 90. The book can be used
as a textbook but it is mainly for reference use. I find it rather
easy to read. It treats and explains everything. It is very similar
to the book of Adams et al.
17. ISO: ISO/IEC 1539:1991, Information Technology - Programming
Languages - Fortran, Second Edition, 1991-07-01, ISO Publications
Department, Case Postale 56, CH-1211 Geneva 20, Switzerland. SFR 185.
The official standard for Fortran 90. Rather difficult as a dictionary.
Requires that you have read a textbook on Fortran 90. The book
is recommended. The standard can also be available in electronic
form both in ASCII and PostScript for a certain charge from Walt
Brainerd, Unicomp Inc., 235 Mt. Hamilton Avenue, Los Altos, CA 94022,
Fax + 1 415 949 4058, E-mail walt@netcom.com.
18. James F. Kerrigan: Migrating to Fortran 90, O'Reilly & Associates,
Sebastopol, CA 1993, 389 pages, $27.95.
Published 19 November 1993. I have not seen it, it claims to be a
practical guide to Fortran 90 for the current Fortran 77 programmer.
19. Elliot B. Koffman and Frank L. Friedman: Problem Solving and
Structured Programming in Fortran 77, Fourth Edition,
Addison-Wesley Publishing Company, Reading, Massachusetts 1990.
ISBN 0-201-51216-5.
A complete textbook in both programming in general and in Fortran
77. Many examples and good exercises. Appendix D treats Fortran
8X (the previous version of Fortran 90). I find this look a
little more easy to read than the one of Ellis. The book is
recommended.
20. Erasmus Langer: Programmieren in Fortran, Springer, Vienna 1993,
ISBN 3-211-82446-4. DEM 45.
Tutorial in German on Fortran 90. Just published, so I have not
read it yet. Contains a unique appendix on the floating point
representation on the most commonly used computers.
21. John M. Levesque and Joel W. Williamson: A Guidebook to Fortran
on Supercomputers, Academic Press, San Diego, CA, 1989. ISBN
0-12-444760-0.
This book treats a lot of tricks in order to vectorize Fortran 77
programs, especially on the Cray. Many of these tricks are however
already included in the Cray compiler. The book also describes some
supercomputer architectures.
22. Mike Loukides: UNIX for Fortran Programmers, Nutshell
Handbooks, O'Reilly & Associates, Sebastopol, CA 1990, ISBN
0-937175-51-X. $24.95.
An excellent UNIX textbook in Fortran programming. It has taught me
how libraries are used in UNIX. The book is recommended.
23. Michael Metcalf: Fortran Optimization, Academic Press, London
and New York 1982. ISBN 0-12-492480-8.
A classical book how you get efficient Fortran 77 programs on a
conventional computer.
24. Michael Metcalf and John Reid: Fortran 90 Explained, Oxford
University Press, Oxford, 1990. ISBN 0-19-853772-7. $29.95.
This book was reprinted with corrections in 1993. A good and rather
easy to read textbook written by persons involved in the
----------------------------------------------------------------------
- 41 -
development of Fortran 90. The 1993 printing contains a very
complete application example.
25. NAG: NAGWare f90 Compiler (Unix), Release 2.0, NP2563, March
1993. ISBN 1-85206-087-5.
A short description of the NAG compiler with the listing of all
Fortran 90 commands and the intrinsic functions. It also contains
some extensions to the standard, three complete modules, and
information on mixing Fortran 90 and C.
26. NAG: FTN90 Reference Manual, NAGWare FTN 90 Professional Plus,
August 1994. ISBN 1-85206-103-0.
A description of NAG's compiler, linker and other utilities.
In addition input/output and modules in Fortran 90 are
discussed. This PC version handbook is much more complete than the
one for UNIX (reference 25).
27. Rama N. Reddy and Carol A. Ziegler: FORTRAN 77 with 90:
Applications for Scientists and Enginners, Second Edition,
West Publishing Company, Minneapolis, 1994. ISBN 0-314-02861-7.
Basically a textbook on Fortran 77 with Fortran 90 extensions
at the end of each chapter.
28. Patrick D. Terry: FORTRAN From Pascal, Addison-Wesley,
Wokingham, England, 1987. ISBN 0-201-17821-4.
The purpose of this book is to be a textbook in Fortran 77 for the
one who knows Pascal. Regrettably, it has more become a book on how
to write such programs, that are in reality more suited for Pascal,
in Fortran 77, e.g. simulation of recursion. Fortran ought to
be used at what it is good for, large numerical or technical
calculations.
29. Christoph Ueberhuber and Peter Meditz: Software-Entwicklung in
Fortran 90, Springer, Vienna 1993, ISBN 3-211-82450-2. DEM 60.
The first part of this book in German discusses the foundations
of numerical computing, and the second part describe Fortran 90.
APPENDICES
APPENDIX 1. Fortran and Pascal
==============================
Advantages of Fortran:
o Fortran is a simple language
o Fortran has always existed
o Fortran compilers are generally available
o Earlier the first programming language
o Used commercially for technical and scientific computations
o Scientists have no interest in learning new languages
o Good at numerical analysis and technical calculations
o It is necessary to structure the problem in order to use Fortran
o A large number of programs and routines in Fortran are exchanged
internationally
o Efficient compilers
o The first standardized programming language
o Better standard obedience than other languages
o Is continually developed (a new version each decade)
o The dominating language on supercomputers
----------------------------------------------------------------------
- 42 -
Differences between Fortran 77 and Pascal:
o Fortran doesn't use assignments with := or end of statements
with ;
o Fortran doesn't have reserved words, it has short identifiers,
the identifiers do not have to be specified
o Fortran does not have records, pointers, user-defined types,
scalar types, subintervals, but it has COMPLEX and DOUBLE
PRECISION
o Fortran 77 does not have WHILE and REPEAT
o Fortran did not get IF THEN ELSE END IF until 1978
o A bad CASE in Fortran 77. Fortran 77 is not able to nest
functions and subroutines and does not permit recursive calls
o Fortran has very good input and output, but those facilities
are very difficult to learn
o Fortran has separate compilation
o Fortran manages national characters in comments and output
o The punched card orientation with Fortran can still give some
problems
o Blanks are not significant (except in the fixed form of
Fortran 90)
o Mixing of integers and floating-point numbers is implemented differently
o Arrays in Fortran 77 have to be assigned values using an explicit loop
o Fortran programs are usually less well structured then Pascal
programs
o Argument association is different
o Common data are treated differently
----------------------------------------------------------------------
- 43 -
APPENDIX 2. Summary of Fortran 77 statements
============================================
Those statements we do not recommend have been indicated with the
question mark "?" and in serious cases even with two question marks
"??".
[ The symbol ? should be replaced by a dagger in the final version ]
Specification of program units:
PROGRAM - main program
FUNCTION - function, FUNCTION can be preceded by some
of the specifications of the variables below,
except IMPLICIT
SUBROUTINE - subroutine
??ENTRY - extra entry in subprograms
? BLOCK DATA - common data, usually given initial values
Specification of variables:
IMPLICIT - default IMPLICIT REAL(A-H, O-Z), INTEGER(I-N)
IMPLICIT NONE - not standard, but very useful, it is available in
Fortran 90. Gives the "Pascal convention" that
all variables have to be specified.
For Sun and DEC the same effect can be obtained
with the switch -u in the compilation command
INTEGER
REAL
DOUBLE PRECISION
COMPLEX
LOGICAL
CHARACTER CHARACTER*4
Additional specifications:
DIMENSION - can also be given directly in the type specification,
as well as in a COMMON
? COMMON - common storage area for variables that are
in several program units
??EQUIVALENCE - common storage area for several variables in the
same program unit
PARAMETER - makes a variable into a constant with a certain value
EXTERNAL - tells the system that the identifier is an
external function or an external subroutine
----------------------------------------------------------------------
- 44 -
INTRINSIC - tells the system that the identifier is an
intrinsic function (or a subroutine, only in
Fortran 90)
SAVE - saves the values between exit or return from one
subroutine into the new call of the same
subroutine or function
DATA - puts initial values into variables
Executable GOTO statements:
GOTO snr1 - ordinary GOTO statement (jumps to the statement with
number snr1)
? GOTO (snr1, snr2, snr3), integer_expression
- conditional GOTO statement. If the integer
expression is 1, 2 or 3, execution jumps to
statement number snr1, snr2 or snr3 (an arbitrary
number of statement numbers snr are permitted).
??GOTO statement_number_variable, (snr1, snr2, snr3)
- an assigned GOTO statement, jumps to the statement
number that equals the statement
number variable (an arbitrary number of
statement numbers snr are permitted).
??GOTO statement_number_variable
- this is an assigned ordinary GOTO statement, it is a
combination of the first one, GOTO snr1, and
previous one, GOTO statement_number_variable without
a list of permitted alternatives.
??ASSIGN statement_number TO statement_number_variable
- statement number variables can not be assigned with
an ordinary assignment of the type (integer
variable = integer expression), it has to be
done with the ASSIGN statement. The statement
number variable can then be used for an assigned
GOTO statement and in the ordinary GOTO statement
and also in connection with FORMAT.
? IF (numerical_expression) snr1, snr2, snr3
- arithmetical IF-statement, jumps to statement number
snr1 if the expression is negative,
snr2 if the expression is zero,
snr3 if the expression is positive
----------------------------------------------------------------------
- 45 -
Other executable statements:
IF(logical_expression) statement
- conditional statement: if the logical expression
is true, the statement is performed, in the
other case execution jumps directly to the next
statement. The statement here is permitted to
be an ordinary assignment statement or an
ordinary jump statement (GOTO statement) or a
call of a subroutine.
IF(logical_expression) THEN ! Complete alternative statement.
...statements... ! Variants without the ELSE-part as well
ELSE ! as with nested ELSE, or with
...statements... ! ELSE replaced by
ENDIF ! ELSE IF (log_expr) THEN
! also exist.
CONTINUE - continuation, does nothing. It is recommended for
clean conclusion of a DO-loop.
STOP - concluding statement, stops execution.
END - concluding statement, stops compilation of the
program unit and also execution if it is
in the main program.
If END is found during execution of a subprogram,
an automatic return to the calling program unit
is executed (replaces the explicit RETURN statement).
? PAUSE - pause statement, stops execution temporarily
(implementation dependent).
DO statement_number variable = var1, var2, var3
- DO-loop.
? Floating-point numbers are permitted as variables
in the DO-loop, but they are not recommended.
It is preferably to use integers.
----------------------------------------------------------------------
- 46 -
Input/output statements:
OPEN - open a file before the program can use it.
CLOSE - close a file. A file that has not been closed can
usually not be read.
READ - input
WRITE - output
PRINT - earlier output to line printer, now a synonym to
WRITE. it works on a standard unit.
INQUIRE - inquires about file status.
REWIND - rewinds a file to the beginning.
BACKSPACE - rewinds a file one record.
ENDFILE - marks end of file.
FORMAT - Fortran speciality (see below).
Call statements:
CALL sbrtn - call a subroutine sbrtn.
fnctn - a function is called by giving the function
name fnctn.
RETURN - return from the subprogram (subroutine or
function).
FORMAT-letters:
Example Comments
-------------------------------------------------------------------
Integer I I5 5 positions reserved
-------------------------------------------------------------------
Floating-point F F8.3 8 positions, out of which 3 are
number used for the fractional part
E E14.6 14 positions of which
6 are used for the decimals
4 - for the exponent
1 - for the sign
1 - for the starting zero
1 - for the decimal point
1 - for a blank character
D D20.12 as E, but for double precision
G G14.6 as F, if the number can be given
within the field, else as E
-------------------------------------------------------------------
Complex numbers as a pair of floating-point variables
-------------------------------------------------------------------
Logical L L1
-------------------------------------------------------------------
----------------------------------------------------------------------
- 47 -
Character A A7 7 characters are available in A7
string
' ' 'Example' Conventional character constant
? nH 7HExample Hollerith constant
Positioning Tn n positions from the left
TLn n positions towards left
TRn n positions towards right
nX n positions towards right
??Last new line $ this is used if you wish to
do input in direct connection
with an output, to stay on
the same line. Not standard!
Not Fortran 90!
Discontinue : if the list does not contain
any more elements the
format is also finished here
New record / normally a new line
-------------------------------------------------------------------
Binary B not Fortran 77 but Fortran 90
Octal O not Fortran 77 but Fortran 90
Hexadecimal Z not Fortran 77 but Fortran 90
-------------------------------------------------------------------
Output SP + is written
SS + is not written
S standard (normal SS)
In all alternatives a minus - is written for negative values
-------------------------------------------------------------------
Input BZ blanks are interpreted as
zeroes
BN blanks are not regarded as
anything (blanks are skipped)
BN is standard using the ULTRIX, when punched cards were used,
BZ was the standard. Compare with BLANK = "ZERO" and
Blank = "NULL" in the OPEN-statement.
-------------------------------------------------------------------
Scaling factor kP:
Input: with an exponent, no action.
Without exponent, the number
is multiplied by 10**(-k)
before assignment, which means
a change of the value.
Output: with exponent, the mantissa
is multiplied by 10**k and the
exponent is reduced with k,
which means no change of the
value.
Without exponent, the
number is multiplied by 10**k
before the output, which means
a change of value.
-------------------------------------------------------------------
NB! S, SP, SS, BN, BZ and kP are valid until the end of the FORMAT or
until a new one of the same kind appears. To scale with kP is
good with E-format on output, because then you avoid that the first
digit is zero, and you get more information into less space on
the paper. To scale with kP is catastrophic using F-format, but
it was of great interest when punched cards were still in use.
A very good and complete description of input and output in Fortran,
including the use of the FORMAT-letters, is given in the book by Adams
et al (1992).
----------------------------------------------------------------------
- 48 -
Extension in Fortran 90 regarding the FORMAT-letters:
-----------------------------------------------------
You can now replace E as the mark for output in exponential form by
ES and then you get the Scientific form with output of one digit
different from zero before the decimal point. If you instead replace E
by EN you get an ENgineering form with one to three digits before the
decimal point and the exponent evenly divisible by three. If the output
value is zero you get the same output from ES and EN as from E.
Addition regarding input/output:
--------------------------------
Most of these have a large number of parameters in the so-called
control list which has been expanded considerably in Fortran 90. They
are treated shortly in five pages of NAG (1992) and are
ACCESS, ACTION, ADVANCE, APPEND, APOSTROPHE, ASIS, BLANK, DELETE,
DELIM, DIRECT, END, EOR, ERR, EXIST, FILE, FMT, FORM, FORMATTED,
IOLENGTH, IOSTAT, KEEP, NAME, NAMED, NEXTREC, NEW, NML, NO, NULL,
NUMBER, OLD, OPENED, PAD, POSITION, QUOTE, READ, READWRITE, REC,
RECL, REPLACE, REWIND, SCRATCH, SEQUENTIAL, SIZE, STATUS,
UNDEFINED, UNFORMATTED, UNIT, UNKNOWN, WRITE, YES and ZERO.
Additions regarding the intrinsic numerical functions:
------------------------------------------------------
The following are available:
ABS, AIMAG, AINT, ANINT, CMPLX, CONJG, DBLE, DIM, DPROD, INT,
MAX, MIN, MOD, NINT, REAL and SIGN.
In addition, CEILING, FLOOR and MODULO have been added to Fortran
90. Only the last one is difficult to explain, which is most easily
done with the examples from ISO (1991)
MOD (8,5) gives 3 MODULO (8,5) gives 3
MOD (-8,5) gives -3 MODULO (-8,5) gives 2
MOD (8,-5) gives 3 MODULO (8,-5) gives -2
MOD (-8,-5) gives -3 MODULO (-8,-5) gives -3
Additions concerning the intrinsic mathematical functions:
----------------------------------------------------------
The following are available:
ACOS, ASIN, ATAN, ATAN2, COS, COSH, EXP, LOG, LOG10, SIN, SINH,
SQRT, TAN and TANH.
Addition concerning output on paper:
------------------------------------
In order to provide a possibility to request for example the top of
the next page there has for many years existed "ASA carriage control
characters", which are described in ISO (1991), section 9.4.5, and in
ANSI (1978), section 12.9.5.2.3. Using UNIX most output is to a
monitor, and the ASA characters are not of interest. If you wish to use
the control characters for output on paper you first have to write on a
file, which is then printed with the tool "fpr" if you are using DEC,
Hewlett Packard or Sun, but a tool called "asa" if you use Cray or
Silicon Graphics. It is used in the following way
----------------------------------------------------------------------
- 49 -
fpr < output_file | lpr -Pprinter
where output_file is the name of the file to be printed with control
characters and printer is the name of the printer to be used.
APPENDIX 3. Summary of the new features in Fortran 90
=====================================================
Below follows a short summary of the new standard. Please note that
only some more essential parts and not the whole new standard are
discussed here.
1. Source code form: As an alternative to the old source-code
(punched card oriented) form, there is a free form, which does not take
any consideration to columns and where blanks are significant. Comments
can be given on the line if preceded by an exclamation mark !, several
statements can be given on the same line if using the semicolon ; as a
separator, and statements that continue on the next line are continued
on the present line (not the next line) with an ampersand &.
The underline symbol _ is permitted inside the names of variables,
and the length of variables must have at most 31 characters (instead of
at most 6 in the Fortran 77 standard).
Blanks become significant in the free form of the source code. Old
commands like ENDIF and GOTO can also in the future be written either as
END IF or GO TO respectively, but of course not like EN DIF or GOT O.
To permit significant blanks in the old (fixed) form of the source code
would not be possible, since it for example is permitted to write END in
the following silly way
E N D
INCLUDE can be used to include source code from an external file.
The construct is a line where there is INCLUDE and a character string
and, perhaps, some concluding comments. Interpretation is
implementation-dependent, normally the character string is treated as
the name of the file that holds the source code that should be included.
Nesting is permitted (and the number of levels is
implementation-dependent), but recursion is not permitted for the
INCLUDE statement.
As in most Algol-type languages the word END can be complemented with
the name of the routine or function like END FUNCTION GAMMA.
2. Alternative representations: The special characters < and > can
now be used
< .LT. > .GT.
<= .LE. >= .GE.
== .EQ. /= .NE.
3. Specifications: These can now be written on one line
REAL, DIMENSION (3), PARAMETER :: &
a = (/ 0.0, 0.0, 0.0 /), b = (/ 1.0, 1.0, 1.0 /)
COMPLEX, DIMENSION(10) :: john
----------------------------------------------------------------------
- 50 -
while the variables a and b become constant vectors with 3 elements
and the floating-point values 0.0 and 1.0, respectively, while john
becomes a complex vector with 10 complex elements, not yet assigned any
values.
If you wish to use the Algol principle to specify all variables, this
is simplified by the command IMPLICIT NONE which switches off the
implicit-type rules.
Double precision has been implemented with a more general method to
give the desired precision, namely the parameter KIND for which
precision we wish, useful on all variable types.
INTEGER, PARAMETER :: LP = SELECTED_REAL_KIND(20)
REAL (KIND = LP) :: X, Y, Z
The above two statements thus declare the variables X, Y and Z to be
REAL floating-point variables with at least 20 decimal digits accuracy
with a data type that is called LP (where LP stands for LONG PRECISION).
4. Conditional statement: The new command is
SELECT CASE (expression)
CASE block-switch
block
CASE block-switch
block
CASE DEFAULT
default block
END SELECT
Typical construct:
SELECT CASE(3*I-J) ! the control variable is 3*i-j
CASE(0) ! for the value zero
: ! you execute the code here
CASE(2,4:7) ! for the variables 2, 4, 5, 6, 7
: ! you execute the code here
CASE DEFAULT ! and for all other values
! you execute the code here
: !
END SELECT
If the CASE DEFAULT is missing and none of the alternatives is valid,
the execution continues directly with the next statement following the
END SELECT, without any error message.
Another example:
INTEGER FUNCTION SIGNUM(N)
SELECT CASE (N)
CASE (:-1)
SIGNUM = -1
CASE (0)
SIGNUM = 0
CASE (1:)
SIGNUM = 1
END SELECT
END
----------------------------------------------------------------------
- 51 -
5. DO-loop: Three new constructs are included. The first one gives
in principle an infinite loop, which however can be terminated with a
conditional GOTO-statement.
name: DO
executable statements
END DO name
The usual DO-loop has the following new simplified form without
statement number,
name: DO i = integer_expr_1, integer_expr_2 ,integer_expr_3
executable statements
END DO name
where i is called control variable, and where ",integer_expr_3" is
optional. Finally there is also the DO WHILE loop
name: DO WHILE (logical_expression)
executable statements
END DO name
The name is optional but can be used for nested loops in order to
indicate which one that is to be iterated once again with the CYCLE
statement or terminated with the EXIT statement.
S1: DO
IF (X > Y ) THEN
Z = X
EXITS1
END IF
CALL NEW(X)
END DO
N = 0
LOOP1: DO I = 1, 10
J= I
LOOP2: DO K =1, 5
L = K
N = N +1
END DO LOOP2
END DO LOOP1
In the latter case the final values from the variables will be as
follows, in full accordance with the standard, I = 11, J = 10, K = 6, L
= 5, and N = 50.
To name the loop is completely optional. Also note that this type of
name is limited to DO-loop, CASE or IF...THEN...ELSE...ENDIF constructs.
The old possibilities with statement numbers are still available, also
in the free form.
6. Program units: Routines can be called with keyword arguments and
can use default arguments
SUBROUTINE solve (a, b, n)
REAL, OPTIONAL, INTENT (IN) :: b
----------------------------------------------------------------------
- 52 -
can be called with
CALL solve (n = i, a = x)
where two of the arguments are given with keywords instead of
position and where the third one has a default value. If SOLVE is an
external routine it requires making use of an INTERFACE block in the
calling program. Routines can be specified to be recursive.
RECURSIVE FUNCTION factorial (n) RESULT (fac)
but must then have a special RESULT name in order to return the
result.
7. Character string variables: CHARACTER has been expanded to
include also an empty string
a = ''
and assignment of an overlapping string is now permitted
a(:5) = a(3:7)
The new intrinsic function TRIM which removes concluding blanks is an
important addition. You can now make a free choice between the
apostrophe ' and the quotation mark " in order to indicate a character
string. This can among other things be used in such a way that if you
wish to write an apostrophe inside the text, then you use the quotation
mark as indicators, and in the opposite case if you wish to have a
quotation mark inside the text you use the apostrophe as the indicator.
8. Input: At last the NAMELIST is included in the standard! This
statement, however, has to be among the specifications. In the example
below list2 is the name of the list, a and b are real variables and i is
an integer variable.
NAMELIST /list2 / a, i, x
:
READ (unit, NML = list2)
which wishes to get input data of the following form, but all
variables do not have to given, and they can be given in any order.
&list2 X = 4.3, A = 1.E20, I = -4 /
9. Vector and matrix management: This is one of the most important
parts of the new standard. An array is defined to have a shape given by
its number of dimensions, called "rank", and the extent for each one of
these. Two arrays agree if they have the same shape. Operations are
normally done element by element. Please remember that the rank for an
array is the number of dimensions and has nothing at all to do with the
mathematical rank of a matrix!
REAL, DIMENSION(5,20) :: x, y
REAL, DIMENSION(-2:2,20) :: z
:
z = 4.0*y*sgrt(x)
We perhaps here wish to protect against negative elements of X. This
is done with the following construct
----------------------------------------------------------------------
- 53 -
WHERE ( x >= 0.0 )
z = 4.0*y*sgrt(x)
ELSEWHERE
z = 0.0
END WHERE
Please note that ELSEWHERE has to be in one word! Compare also with
the function SUM which is discussed at the end of the next section.
You can pick out a part of an array. Assume that the array A is
specified in the following way.
REAL, DIMENSION(-4:0, 7) :: A
With A(-3, :) you pick the second row, while with A(0:-4:-2, 1:7:2)
you pick (in reverse order) its each other element in each other column.
Just as variables can form arrays, also constants can form arrays.
REAL, DIMENSION(6) :: B
REAL, DIMENSION(2,3) :: C
B = (/ 1, 1, 2, 3, 5, 8 /)
C = RESHAPE( B, (/ 2,3 /) )
where the first argument to the intrinsic function RESHAPE gives the
value and the second argument gives the new shape. Two additional, but
optional, arguments are available to this function.
The above can also be written in a more compressed form using the
PARAMETER attribute. In the first line below the PARAMETER attribute is
compulsory (if the assignment is to be made on the same line), but in
the second line it is optional. Remember that the PARAMETER attribute
means that the quantity can not be changed during execution of the
program.
REAL, DIMENSION(6), PARAMETER :: B = (/ 11, 12, 13, 14, 15, 16 /)
REAL, DIMENSION(2,3), PARAMETER :: C = RESHAPE( B, (/ 2, 3 /) )
Any statements for real parallel computation are not included in
Fortran 90. The committee believes it is necessary with additional
experience before the standardization of parallelization. See also HPF
discussed in the Appendix 10.
10. Dynamic storage: Fortran 90 contains four different ways to make
dynamical access. The first one is to use a pointer. See exercise
(12.3)
The second is to use an "allocatable array", i.e. with the
statements ALLOCATE and DEALLOCATE you get and return a storage area for
an array with type, rank and name (and possible other attributes) which
had been specified earlier with the additional attribute ALLOCATABLE.
REAL, DIMENSION(:), ALLOCATABLE :: x
:
Allocate(x(N:M)) ! N and M are the integer expressions here.
:
x(j) = q ! Some assignment of the array.
CALL sub(x) ! Use of the array in a subroutine.
:
DEALLOCATE (x)
----------------------------------------------------------------------
- 54 -
Deallocation occurs automatically (if the attribute SAVE has not been
given) when you reach RETURN or END in the same program unit.
The third variant is an "automatic array", it is almost available in
the old Fortran, where x in the example below has to be in the list of
arguments. This is not required any more.
SUBROUTINE sub (i, j, k)
REAL, DIMENSION (i, j, k) :: x
Dimensions for x are taken from the integers in the calling program.
Finally there is an "assumed-shape array" where the storage is
defined in the calling procedure and for which only the type, rank and
name are given.
SUBROUTINE sub(a)
REAL, DIMENSION (:,:,:) :: a
According to Metcalf and Reid (1990, 1992), section 6.3 you here
require an explicit interface. This has to look as follows
INTERFACE
SUBROUTINE SUB(A)
REAL, DIMENSION (:,:,:) :: A
END SUBROUTINE SUB
END INTERFACE
If you forget the INTERFACE or if you have an erroneous interface,
then you will usually get "segmentation error", it means that a program
unit may be missing. See also the discussion on page 15.
Some intrinsic functions are available to determine the actual
dimension limits
DO (i = LBOUND(a,1), UBOUND(a,1))
DO (j = LBOUND (a,2), UBOUND (a,2))
DO (k = LBOUND(a,3),UBOUND (a,3))
where LBOUND gives the lower limit for the specified dimension and
UBOUND gives the upper one.
The sum of the positive value of a number of elements in an array is
written
SUM ( X, MASK = X .GT. 0.0)
These statements can not be used in order to avoid division by zero
at for example summation of 1/X, that is the mask works only with
determining which numbers that are to be included in the summation, and
not whether a certain value has to be calculated or not. But in this
later case you can use the construct WHERE, see section 9.
11. Intrinsic functions: Fortran 90 defines about 100 intrinsic
functions and a few intrinsic subroutines. Many of these functions can
be used for arrays, e.g. for reduction (SUM), construct (SPREAD),
manipulation (TRANSPOSE). Other permits attributes or available
parameters in the programming environment to be determined, as the
largest positive floating-point number and the largest positive integer,
as well as access to the system clock. The random numbers generator is
included. Finally, the function TRANSFER permits a certain physical
area to be transmitted to another area without type conversion.
----------------------------------------------------------------------
- 55 -
The function SPREAD is discussed more fully in the solution of
exercise (11.1).
All intrinsic functions and subroutines are discussed in Appendix 5.
12. Use of the user-defined data type: Fortran had earlier not
permitted use of any user-defined type. This has now become possible.
TYPE staff_member
CHARACTER(LEN=20) :: first_name, last_name
INTEGER :: identification, department
END TYPE
which can be used in order to describe an individual. A combination
of individuals can also be formed
TYPE(staff_member), DIMENSION(100) :: staff
Individuals can be referred to as staff(number) and a field can be
referred as staff(number)%first_name. You can also nest definitions
TYPE company
CHARACTER(LEN=20) :: company_name
TYPE(staff_member), DIMENSION(100) :: staff
END TYPE
:
TYPE(company), DIMENSION(10) :: several_companies
A numerically more interesting example is a sparse matrix A with at
most one hundred non-zero elements, which can be specified with the
following statement
TYPE NONZERO
REAL VALUE
INTEGER ROW, COLUMN
END TYPE
and
TYPE (NONZERO) :: A(100)
You then get the value of A(10) by writing A(10)%VALUE. Assignment
can be done, for example with A(15) = NONZERO(17.0,3,7).
In order to use user-defined data types in for example COMMON, or to
make sure that two data types which look the same are treated as
identical, you can use the SEQUENCE statement, in the latter case it is
also required that no variable is specified PRIVATE.
13. Modules are collections of data, type definitions and procedure
definitions, which give a more secure and general replacement for the
COMMON concept.
14. The data type "bit" is not included in the standard, but there
are available various bit operations of integers according to an earlier
military standard MIL-STD 1753. In addition, binary, octal and
hexadecimal constants are included, as well as the possibility to use
these quantities in input/output through the three new format-letters.
In a DATA statement you can use an assignment to
B'01010101010101010101010101010101'
----------------------------------------------------------------------
- 56 -
for binary,
O'01234567'
for octal, and
Z'ABCDEF'
for hexadecimal numbers.
15. Pointers have been included, but not in the usual way as a new
data type but as an attribute to the other data types. A variable with
a pointer attribute can be used as an ordinary variable and in a number
of new ways. Pointers in Fortran 90 are not memory addresses as in many
languages or in certain Fortran variants (dialects), but are more like
extra names (aliases).
16. User extensions: The new language contains a possibility for the
user to extend it with his own concepts, for example interval
arithmetics, rational arithmetics or dynamic character strings. By
defining a new data type or operator, and overloading operations and
procedures (so that you can also use the plus + as the symbol of
addition of intervals and not only of ordinary numbers), we can create a
package (a module) without using a preprocessor. We can soon expect a
number of extensions for different applications in the form of modules
from different manufacturers. Some are already available from NAG.
----------------------------------------------------------------------
- 57 -
APPENDIX 4. Backward and Forward compatibility
==============================================
Very important in the introduction of a new programming language
standard is that old programs (at least those who obey the outgoing
standard) can be used once again, with the new standard.
When we went from Fortran 66 to Fortran 77 the extended DO-loop was
removed (the extended DO-loop means that if you do not change any of the
DO-loop parameters you can jump out of the loop and then jump in again
(this is somewhat contrary to the concept of structured programming).
In addition Hollerith constants were removed (except in FORMAT). That
means that there are some programs that obey Fortran 66 but do not obey
Fortran 77. Most manufacturers have, however, chosen to let these two
concepts be included as extensions in their Fortran implementations.
For Fortran 90 nothing has been removed from Fortran 77. An interesting
practical question is however manufacturers still continue to include
those old things that really should have been thrown away when Fortran
77 came. It is permitted to have these old concepts as extensions.
On the other hand, the concept "obsolescence" is introduced. This
means that some constructs may be removed at the next change of Fortran.
These constructs are:
Arithmetic IF-statement
Control variables in a DO-loop which are floating point or
double-precision floating-point
Terminating several DO-loops on the same statement
Terminating the DO-loop in some other way than with CONTINUE
or END DO
Alternate return
Jump to END IF from an outer block
PAUSE
ASSIGN and "assigned GOTO" and "assigned FORMAT", that is the whole
"statement number variable" concept.
Hollerith editing in FORMAT
A group "High Performance Fortran Forum" has worked at the
development of an extension to Fortran 90 with parallel extensions. The
purpose of this project is to offer a portable language which can be
used efficiently on different parallel systems. The project was ready,
with a complete proposal, in May 1993 and aims at a de facto standard
(and not at a formal standard). See Appendix 10 for a summary.
Somewhat simplified you can say that Fortran 90 works effectively on
vector processors but not on parallel processors.
Among the new things that are being considered for the next version
of Fortran are improved parallel treatment, interrupt handling,
parametrized data types, and data types with inherited properties. It
is the aim of the committee to have a slight revision in 1996, with some
carefully chosen new properties.
In addition, a few corrections will come sooner, especially some
explanations of ambitious parts of the standard. An early version of
these corrections appeared in a special issue of the Fortran Forum, Vol.
11, No. 1, March 1993, with the title "Fortran 90; Errata, Amendments,
and Interpretations, progress to date". Some updates to this special
issue appeared in No. 2, June 1993, page 1.
----------------------------------------------------------------------
- 58 -
APPENDIX 5. Intrinsic functions in Fortran 90
=============================================
There is a large a number of intrinsic functions and five intrinsic
subroutines in Fortran 90. I treat the numeric and mathematical
routines very shortly, since they are not changed from Fortran 77 and
therefore should be well-known.
This section is based on section 13 of the ISO standard (1991), which
contains a more formal treatment. We follow the arrangement of the
different functions and subroutines in the standard, but explain
directly in the list. For a more detailed treatment we refer to Metcalf
and Reid (1990, 1993).
When a parameter below is optional it is given in lower case
characters. When an argument list contains several arguments the
function can be called either by position related arguments or by a
keyword. Keyword must be used if some previous argument is not
included. Keywords are normally the names that are given below.
We have not always given all the natural limitations to the
variables, for example that the rank is not permitted to be negative.
1. Function which determines if a certain argument is in an actual
argument list:
The function
PRESENT(A)
returns .TRUE. if the argument A is in the calling list, .FALSE. in
the other case. The use is illustrated in section 8 of the main text.
2. Numerical functions:
The following functions from Fortran 77 can use a kind-parameter like
in AINT(A, kind), namely AINT, ANINT, CMPLX, INT, NINT and REAL. In
addition, three new functions were introduced with Fortran 90.
See page 48 for a list of all the functions.
[This list will me moved here]
3. Mathematical functions:
Same as in Fortran 77. See page 48 for a list of all the functions.
[This list will me moved here]
4. Character string functions:
The functions below perform operations from and to character strings.
Please note that ACHAR works with the standard ASCII character set while
CHAR works with the representation in the computer you are using.
ACHAR(I) Returns the ASCII character which has number I
ADJUSTL(STRING) Adjusts to the left
ADJUSTR(STRING) Adjusts to the right
CHAR(I, kind) Returns the character that has the number I
IACHAR(C) Returns the ASCII number of the character C
ICHAR(C) Returns the number of character C
----------------------------------------------------------------------
- 59 -
INDEX(STRING, SUBSTRING, back) Returns the starting position for a
substring within a string. If
BACK is true then you get the last
starting position, in the other
case, the first one.
LEN_TRIM(STRING) Returns the length of the string
without the possibly trailing blanks.
LGE(STRING_A, STRING_B)
LGT(STRING-A, STRING_B)
LLE(STRING_A, STRING_B)
LLT(STRING_A, STRING_B)
The above routines compare two strings using sorting according to
ASCII. If a string is shorter than the other, blanks are added at the
end of the short string. If a string contains a character outside the
ASCII character set, the result is implementation-dependent.
[Bo, how do you deal with Swedish characters then ?]
[Reply: With the seven-bit ASCII-code the Swedish characters replace
another character, ae is equivalent to {, aa to } and oe to |,
and AE to [, AA to ] and finally OE to \ . With the eight-bit ASCII
we do not manage, we have to use the implementation dependent functions.]
REPEAT(STRING, NCOPIES) Concatenates a character string NCOPIES
times with itself.
SCAN(STRING, SET, back) Returns the position of the first occurrence
of any character in the string SET in the string
STRING. If BACK is true, you will get
the rightmost such character.
TRIM(STRING) Returns the character string STRING without
trailing blanks.
VERIFY(STRING, SET, back) Returns the position of the first character
in STRING which is not in SET. If BACK
is TRUE, you get the last one!
The result is zero if all characters are
included!
5. Character string function for request:
LEN(STRING) This routine returns the length of a character
string. There does not have to be assigned a
value to the variable STRING.
6. Kind functions:
KIND(X)
SELECTED_INT_KIND(R)
SELECTED_REAL_KIND(p, r)
The first returns the kind of the actual argument, which can be of
the type INTEGER, REAL, COMPLEX, LOGICAL or CHARACTER. The argument X
does not have to be assigned any value. The second returns an integer
kind with the requested number of digits, and the third returns the kind
for floating-point numbers
----------------------------------------------------------------------
- 60 -
with numerical precision at least P digits and one decimal exponent
range between -R and +R. The parameters P and R must be scalar
integers. At least one of P and R must be given.
[which case should be p and r - lower or upper?]
[Reply: I am following the notation that optional parameters are given
with lower case at the definition, but upper case at the actual usage]
The result of SELECTED_INT_KIND is an integer from zero and upward,
if the desired kind is not available you will get -1. If several
implemented types satisfy the condition, the one with the least decimal
range is used. If there still are several types or kinds that satisfy
the condition, the one with the smallest kind number will be used.
The result of SELECTED_REAL_KIND is also an integer from zero and
upward; if the desired kind is not available, then -1 is returned if the
precision is not available, -2 if the exponent range is not available
and -3 if no one of the requirements are available. If several
implemented types satisfy the condition, the one with the least decimal
precision is returned, and if there are several of them, the one with
the least kind number is returned.
Examples are given in Section 2 of the main text. Examples of kinds
in a few different implementations are given in Appendix 6.
7. Logical function:
LOGICAL(L, kind)
Converts between different kinds of logical variables. Logical
variables can be implemented in various ways, for example with a
physical representation occupying one bit (not recommended), one byte,
one word or perhaps even one double word. This difference is important
if COMMON and EQUIVALENCE with logical variables have been misused in a
program in the traditional way of Fortran 66 programming.
8. Numerical inquiry functions:
These functions work with a certain model of integer and
floating-point arithmetics, see ISO (1991), section 13.7.1. The
functions return properties of numbers of the same kind as the variable
X, which can be real and in some cases integer.
Functions that return properties of the actual argument X are
available under item 12 - floating-point manipulation functions
DIGITS(X) The number of significant digits
EPSILON(X) The least positive number that added to 1
returns a number that is greater than 1
HUGE(X) The largest positive number
MAXEXPONENT(X) The largest exponent
MINEXPONENT The smallest exponent
PRECISION(X) The decimal precision
RADIX(X) The base in the model
RANGE(X) The decimal exponent
TINY(X) The smallest positive number
----------------------------------------------------------------------
- 61 -
9. Bit inquiry function:
BIT_SIZE(I) Returns the number of bits according to the
model of bit representation in the standard
ISO (1991), section 13.5.7. Normally we get the
number of bits in a (whole) word.
10. Bit manipulation functions:
The model for bit representation in the standard ISO( 1991), section
13.5.7 is used.
BTEST(I, POS) TRUE if the position number POS of I is 1
IAND(I, J) logical addition of the bit characters in
variables I and J
IBCLR(I, POS) puts a zero in the bit in position POS
IBITS(I, POS, LEN) uses LEN bits of the word I with
beginning in position POS, the additional bits
are set to zero. It requires that
POS + LEN <= BIT_SIZE(I)
IBSET(I, POS) puts the bit in position POS to 1
IEOR(I, J) performs logical exclusive OR
IOR(I, J) performs logical OR
ISHIFT(I, SHIFT) performs logical shift (to the right if the number
of steps SHIFT < 0 and to the left if SHIFT > 0).
Positions that are vacated are set to zero.
ISHIFTC(I, SHIFT, size) performs logical shift a number of steps
circularly to the right if SHIFT < 0,
circularly to the left if SHIFT > 0. If SIZE
is given, it is required that 0 < SIZE <=
BIT_SIZE(I). Shift is only done for the bits
that are in the SIZE rightmost positions, but
for all positions if SIZE is not given.
NOT(I) returns a logical complement
11. Transfer functions:
TRANSFER(SOURCE, MOULD, size)
This function specifies that the physical representation of the first
argument SOURCE shall be treated as if it had type and parameters as the
second argument MOULD, but without converting it. The purpose is to
give a possibility to move a quantity of a certain type via a routine
that does not have exactly that data type.
12. Floating-point manipulation functions:
These functions work in a certain model of integer and floating-point
arithmetic, see the standard ISO(1991), section 13.7.1. The functions
return numbers related to the actual variable X of the type REAL.
Functions that return properties for the numbers of the same kind as the
variable X are under item 8 (Numerical inquiry functions).
EXPONENT(X) exponent of the number
FRACTION(X) the fractional part of the number
NEAREST(X, S) returns the next representable number in the
direction of the sign of S
RRSPACING(X) returns the inverted value of the distance between
the two nearest possible numbers
SCALE(X, I) multiplies X by the base to the power I
----------------------------------------------------------------------
- 62 -
SET_EXPONENT(X, I) returns the number that has the fractional part
of X and the exponent I
SPACING(X) the distance between the two nearest possible
numbers
13. Vector- and matrix-multiplication functions:
DOT_PRODUCT(VECTOR_A, VECTOR_B)
Makes a scalar product of two vectors, which must have the same
length (same number of elements).
MATMUL(MATRIX_A, MATRIX_B)
Makes the matrix product of two matrices, which must be consistent,
i.e. have the dimensions like (M, K) and (K, N). Used in section 11 of
the main text.
14. Array functions:
ALL(MASK, dim)
Returns a logical value that indicates whether all relations in MASK
are TRUE, along only the desired dimension if the second argument is
given.
ANY(MASK, dim)
Returns a logical value that indicates whether any relation in MASK
is TRUE, along only the desired dimension if the second argument is
given.
COUNT(MASK, dim)
Returns a numerical value that is the number of relations in MASK who
are TRUE, along only the desired dimension if the second argument is
given.
MAXVAL(ARRAY, dim, mask)
Returns the largest value in the array ARRAY, of those that obey the
relation in the third argument MASK if that one is given, along only the
desired dimension if the second argument DIM is given.
MINVAL(ARRAY, dim, mask)
Returns the smallest value in the array ARRAY, of those that obey the
relation in the third argument MASK if that one is given, along only the
desired dimension if the second argument DIM is given.p
PRODUCT(ARRAY, dim, mask)
Returns the product of all the elements in the array ARRAY, of those
that obey the relation in the third argument MASK if that one is given,
along only the desired dimension if the second argument DIM is given.
SUM (ARRAY, dim, mask)
Returns the sum of all the elements in the array ARRAY, of those that
obey the relation in the third argument MASK if that one is given, along
only the desired dimension if the second argument DIM is given. An
example is given in Appendix 3, section 10.
----------------------------------------------------------------------
- 63 -
15. Array inquiry functions: See also Appendix 3, item 10.
ALLOCATED(ARRAY)
Logical function which indicates if the array is allocated.
LBOUND(ARRAY, dim)
Function which returns the lower dimension limit for the ARRAY. If
DIM (dimension) is not given as an argument, you get an integer vector,
if DIM is included, you get the integer value with exactly that lower
dimension limit, for which you asked.
SHAPE(SOURCE)
Function which returns the shape of an array SOURCE as an integer
vector.
SIZE(ARRAY, dim)
Function which returns the number of elements in an array ARRAY, if
DIM is not given, and the number of elements in the relevant dimension
if DIM is included.
UBOUND(ARRAY, dim)
Function similar to LBOUND which returns the upper dimensional
limits.
16. Array construct functions:
MERGE(TSOURCE, FSOURCE, MASK)
Function which joins two arrays. It gives the elements in TSOURCE if
the condition in MASK is TRUE and FSOURCE if the condition in MASK is
FALSE. The two fields TSOURCE and FSOURCE have to be of the same type
and the same shape. The result is also of this type and this shape.
Also MASK must have the same shape.
I here give a rather complete example of the use of MERGE which also
uses RESHAPE from the next section in order to build suitable test
matrices.
Note that the two subroutines WRITE_ARRAY and WRITE_L_ARRAY are test
routines to write matrices which in the first case are of a REAL type,
in the second case of a LOGICAL type.
IMPLICIT NONE
INTERFACE
SUBROUTINE WRITE_ARRAY (A)
REAL :: A(:,:)
END SUBROUTINE WRITE_ARRAY
SUBROUTINE WRITE_L_ARRAY (A)
LOGICAL :: A(:,:)
END SUBROUTINE WRITE_L_ARRAY
END INTERFACE
----------------------------------------------------------------------
- 64 -
REAL, DIMENSION(2,3) :: TSOURCE, FSOURCE, RESULT
LOGICAL, DIMENSION(2,3) :: MASK
TSOURCE = RESHAPE( (/ 11, 21, 12, 22, 13, 23 /), &
(/ 2, 3 /) )
FSOURCE = RESHAPE( (/ -11, -21, -12, -22, -13, -23 /), &
(/ 2,3 /) )
MASK = RESHAPE( (/ .TRUE., .FALSE., .FALSE., .TRUE.,&
.FALSE., .FALSE. /), (/ 2,3 /) )
RESULT = MERGE(TSOURCE, FSOURCE, MASK)
CALL WRITE_ARRAY(TSOURCE)
CALL WRITE_ARRAY(FSOURCE)
CALL WRITE_ARRAY(MASK)
CALL WRITE_ARRAY(RESULT)
END
SUBROUTINE WRITE_ARRAY (A)
REAL :: A(:,:)
DO I = LBOUND(A,1), UBOUND(A,1)
WRITE(*,*) (A(I, J), J = LBOUND(A,2), UBOUND(A,2) )
END DO
RETURN
END SUBROUTINE WRITE_ARRAY
SUBROUTINE WRITE_L_ARRAY (A)
LOGICAL :: A(:,:)
DO I = LBOUND(A,1), UBOUND(A,1)
WRITE(*,"(8L12)") (A(I, J), J= LBOUND(A,2), UBOUND(A,2))
END DO
RETURN
END SUBROUTINE WRITE_L_ARRAY
The following output is obtained
11.0000000 12.0000000 13.0000000
21.0000000 22.0000000 23.0000000
-11.0000000 -12.0000000 -13.0000000
-21.0000000 -22.0000000 -23.0000000
T F F
F T F
11.0000000 -12.0000000 -13.0000000
-21.0000000 22.0000000 -23.0000000
----------------------------------------------------------------------
- 65 -
PACK(ARRAY, MASK, vector)
Packs an array to a vector with the control of MASK. The shape of
the logical array MASK has to agree with the one for ARRAY or MASK must
be a scalar. If VECTOR is included, it has to be an array of rank 1
(i.e. a vector) with at least as many elements as those that are true
in MASK and have the same type as ARRAY. If MASK is a scalar with the
value TRUE then VECTOR instead must have the same number of elements as
ARRAY.
The result is a vector with as many elements as those in ARRAY that
obey the conditions if VECTOR is not included (i.e. all elements if
MASK is a scalar with value TRUE). In the other case the number of
elements of the result will be as many as in VECTOR. The values will be
the approved ones, i.e. the values which fulfil the condition, and will
be in the ordinary Fortran order. If VECTOR is included and the number
of its elements exceeds the number of approved values, the lacking
values required for the result are taken from the corresponding
locations in VECTOR.
The following example is based on the modification of the one for
MERGE, but I give now only the results.
ARRAY
11.0000000 12.0000000 13.0000000
21.0000000 22.0000000 23.0000000
VECTOR
-11.0000000
-21.0000000
-12.0000000
-22.0000000
-13.0000000
-23.0000000
MASK
T F F
F T F
PACK(ARRAY, MASK)
11.0000000
22.0000000
PACK(ARRAY, MASK, VECTOR)
11.0000000
22.0000000
-12.0000000
-22.0000000
-13.0000000
-23.0000000
----------------------------------------------------------------------
- 66 -
SPREAD(SOURCE, DIM, NCOPIES)
The function SPREAD (SOURCE, DIM, NCOPIES) returns an array of the
same type as the argument SOURCE with the rank increased by one. The
parameters DIM and NCOPIES are integer. If NCOPIES is negative the
value zero is used instead. If SOURCE is a scalar, then SPREAD becomes
a vector with NCOPIES elements that all have the same value as SOURCE.
The parameter DIM indicates which index is to be extended. It has to be
within the range 1 and 1+(rank of SOURCE), if SOURCE is a scalar then
DIM has to be one. The parameter NCOPIES is the number of elements in
the new dimensions. Additional discussion is given in the solution to
exercise (11.1).
UNPACK(VECTOR, MASK, ARRAY)
This function scatters a vector to an array under control of MASK.
The shape of the logical array MASK has to agree with the one for ARRAY.
The array VECTOR has to have the rank 1 (i.e. it is a vector) with at
least as many elements as those that are true in MASK, and also has to
have the same type as ARRAY. If ARRAY is given as a scalar then it is
considered to be an array with the same shape as MASK and the same
scalar elements everywhere.
The result will be an array with the same shape as MASK and the same
type as VECTOR. The values will be those from VECTOR that are accepted
(i.e. those fulfilling the condition in MASK), taken in the ordinary
Fortran order, while in the remaining positions in ARRAY the old values
are kept.
17. ARRAY reshape function.
RESHAPE(SOURCE, SHAPE, pad, order)
Constructs an array with a specified shape SHAPE starting from the
elements in a given array SOURCE. If PAD is not included then the size
of SOURCE has to be at least PRODUCT (SHAPE). If PAD is included it has
to have the same type as SOURCE. If ORDER is included, it has to be an
INTEGER array with the same shape as SHAPE and the values must be a
permutation of (1,2,3,...,N), where N is the number of elements in
SHAPE, it has to be less than, or equal to 7.
The result has of course a shape SHAPE and the elements are those in
SOURCE, possibly complemented with PAD. The different dimensions have
been permuted at the assignment of the elements if ORDER was included,
but without influencing the shape of the result.
A few simple examples are given in the previous and the next section
and also in Appendix 3, item 9. A more complicated example,
illustrating also the optional arguments, follows.
! PROGRAM TO TEST THE OPTIONAL ARGUMENTS TO RESHAPE
INTERFACE
SUBROUTINE WRITE_MATRIX(A)
REAL, DIMENSION(:,:) :: A
END SUBROUTINE WRITE_MATRIX
END INTERFACE
REAL, DIMENSION (1:9) :: B = (/ 11, 12, 13, 14, 15, 16, 17, 18, 19 /)
REAL, DIMENSION (1:3, 1:3) :: C, D, E
REAL, DIMENSION (1:4, 1:4) :: F, G, H
----------------------------------------------------------------------
- 67 -
INTEGER, DIMENSION (1:2) :: ORDER1 = (/ 1, 2 /)
INTEGER, DIMENSION (1:2) :: ORDER2 = (/ 2, 1 /)
REAL, DIMENSION (1:16) :: PAD1 = (/ -1, -2, -3, -4, -5, -6, -7, -8, &
& -9, -10, -11, -12, -13, -14, -15, -16 /)
C = RESHAPE( B, (/ 3, 3 /) )
CALL WRITE_MATRIX(C)
D = RESHAPE( B, (/ 3, 3 /), ORDER = ORDER1)
CALL WRITE_MATRIX(D)
E = RESHAPE( B, (/ 3, 3 /), ORDER = ORDER2)
CALL WRITE_MATRIX(E)
F = RESHAPE( B, (/ 4, 4 /), PAD = PAD1)
CALL WRITE_MATRIX(F)
G = RESHAPE( B, (/ 4, 4 /), PAD = PAD1, ORDER = ORDER1)
CALL WRITE_MATRIX(G)
H = RESHAPE( B, (/ 4, 4 /), PAD = PAD1, ORDER = ORDER2)
CALL WRITE_MATRIX(H)
END
SUBROUTINE WRITE_MATRIX(A)
REAL, DIMENSION(:,:) :: A
WRITE(*,*)
DO I = LBOUND(A,1), UBOUND(A,1)
WRITE(*,*) (A(I,J), J = LBOUND(A,2), UBOUND(A,2))
END DO
END SUBROUTINE WRITE_MATRIX
The output from the above program is as follows.
11.0000000 14.0000000 17.0000000
12.0000000 15.0000000 18.0000000
13.0000000 16.0000000 19.0000000
11.0000000 14.0000000 17.0000000
12.0000000 15.0000000 18.0000000
13.0000000 16.0000000 19.0000000
11.0000000 12.0000000 13.0000000
14.0000000 15.0000000 16.0000000
17.0000000 18.0000000 19.0000000
11.0000000 15.0000000 19.0000000 -4.0000000
12.0000000 16.0000000 -1.0000000 -5.0000000
13.0000000 17.0000000 -2.0000000 -6.0000000
14.0000000 18.0000000 -3.0000000 -7.0000000
11.0000000 15.0000000 19.0000000 -4.0000000
12.0000000 16.0000000 -1.0000000 -5.0000000
13.0000000 17.0000000 -2.0000000 -6.0000000
14.0000000 18.0000000 -3.0000000 -7.0000000
11.0000000 12.0000000 13.0000000 14.0000000
15.0000000 16.0000000 17.0000000 18.0000000
19.0000000 -1.0000000 -2.0000000 -3.0000000
-4.0000000 -5.0000000 -6.0000000 -7.0000000
----------------------------------------------------------------------
- 68 -
18. ARRAY manipulation functions.
The shift functions return the shape of an array unchanged, but move
the elements. They are rather difficult to explain so I recommend to
study also the standard ISO (1991).
CSHIFT(ARRAY, SHIFT, dim)
Performs circular shift by SHIFT positions to the left if SHIFT is
positive and to the right if it is negative. If ARRAY is a vector the
shift is being done in a natural way, if it is an array of a higher rank
then the shift is in all sections along the dimension DIM. If DIM is
missing
it is considered to be 1, in other cases it has to be a scalar
integer number between 1 and n (where n equals the rank of ARRAY). The
argument SHIFT is a scalar integer or an integer array of rank n-1 and
the same shape as the ARRAY, except along the dimension DIM (which is
removed because of the lower rank). Different sections can therefore be
shifted in various directions and with various numbers of positions.
EOSHIFT(ARRAY, SHIFT, boundary, dim)
Performs shift to the left if SHIFT is positive and to the right if
it is negative. Instead of the elements shifted out new elements are
taken from BOUNDARY. If ARRAY is a vector the shift is being done in a
natural way, if it is an array of a higher rank, the shift on all
sections is along the dimension DIM. If DIM is missing, it is
considered to be 1, in other cases it has to have a scalar integer value
between 1 and n (where n equals the rank of ARRAY). The argument SHIFT
is a scalar integer if ARRAY has rank 1, in the other case it can be a
scalar integer or an integer array of rank n-1 and with the same shape
as the array ARRAY except along the dimension DIM (which is removed
because of the lower rank). The corresponding applies to BOUNDARY which
has to have the same type as the ARRAY. If the parameter BOUNDARY is
missing you have the choice of values zero, .FALSE. or blank being
used, depending on the data type. Different sections can thus be
shifted in various directions and with various numbers of positions.
A simple example of the above two functions for the vector case
follows, both the program and the output.
REAL, DIMENSION(1:6) :: A = (/ 11.0, 12.0, 13.0, 14.0, &
15.0, 16.0 /)
REAL, DIMENSION(1:6) :: X, Y
WRITE(*,10) A
X = CSHIFT ( A, SHIFT = 2)
WRITE(*,10) X
Y = CSHIFT (A, SHIFT = -2)
WRITE(*,10) Y
X = EOSHIFT ( A, SHIFT = 2)
WRITE(*,10) X
Y = EOSHIFT ( A, SHIFT = -2)
WRITE(*,10) Y
10 FORMAT(1X,6F6.1)
END
11.0 12.0 13.0 14.0 15.0 16.0
13.0 14.0 15.0 16.0 11.0 12.0
15.0 16.0 11.0 12.0 13.0 14.0
13.0 14.0 15.0 16.0 0.0 0.0
0.0 0.0 11.0 12.0 13.0 14.0
----------------------------------------------------------------------
- 69 -
A simple example of the above two functions in the matrix case
follows. I have here used RESHAPE in order to create a suitable
matrix to start work with. The program is not reproduced here, only
the main statements.
B = (/ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 /)
11.0 12.0 13.0 Z = RESHAPE( B, (/3,3/) )
14.0 15.0 16.0
17.0 18.0 19.0
17.0 18.0 19.0 X = CSHIFT (Z, SHIFT = 2)
11.0 12.0 13.0
14.0 15.0 16.0
13.0 11.0 12.0 X = CSHIFT ( Z, SHIFT = 2, DIM = 2)
16.0 14.0 15.0
19.0 17.0 18.0
14.0 15.0 16.0 X = CSHIFT (Z, SHIFT = -2)
17.0 18.0 19.0
11.0 12.0 13.0
17.0 18.0 19.0 X = EOSHIFT ( Z, SHIFT = 2)
0.0 0.0 0.0
0.0 0.0 0.0
13.0 0.0 0.0 X = EOSHIFT ( Z, SHIFT = 2, DIM = 2)
16.0 0.0 0.0
19.0 0.0 0.0
0.0 0.0 0.0 X = EOSHIFT ( Z, SHIFT = -2)
0.0 0.0 0.0
11.0 12.0 13.0
TRANSPOSE (MATRIX)
Transposes a matrix, which is an array of rank 2. It replaces the
rows and columns in the matrix.
19. Array location functions:
MAXLOC(ARRAY, mask)
Returns the position of the greatest element in the array ARRAY, if
MASK is included only for those which fulfil the conditions in MASK.
The result is an integer vector! It is used in the solution of exercise
(11.1).
MINLOC(ARRAY, mask)
Returns the position of the smallest element in the array ARRAY, if
MASK is included only for those which fulfil the conditions in MASK.
The result is an integer vector!
----------------------------------------------------------------------
- 70 -
20. Pointer inquiry functions:
ASSOCIATED(POINTER, target)
A logical function that indicates if the pointer POINTER is
associated with some target, and if a specific TARGET is included it
indicates if it is associated with exactly that target. If both POINTER
and TARGET are pointers, the result is TRUE only if both are associated
with the same target.
I refer the reader to section 12 of the main text, Pointers.
21. Intrinsic subroutines:
Time routines:
DATE_AND_TIME(date, time, zone, values)
A subroutine which returns the date, the time and the time zone. At
least one argument has to be given.
DATE must be a scalar character string variable with at least 8
characters and it is assigned the value CCYYMMDD for century, year,
month and day. All are given numerically, with blanks if the system
does not include the date.
TIME must also be a scalar character string variable with at least 10
characters and it is assigned a value hhmmss.sss for time in hours,
minutes, seconds and milliseconds. All are given numerically with
blanks if the system does not include a clock.
ZONE must be a scalar character string variable with at least 5
characters and it is assigned the value +hhmm for sign, time in hours
and minutes for the local time difference with UTC (which was previously
called Greenwich Mean Time). All are given numerically, with blanks if
the system does not include a clock. In Sweden we therefore get +0100
in winter and +0200 in summer, in Novosibirsk we get +0700.
[Is the time difference towards GMT the same in the summer and winter
in Russia, i.e. do you have summer time or Daylight Saving Time? / Bo]
The variable VALUES is instead an integer vector with at least 5
elements, it gives the easiest way of using the results from
DATE_AND_TIME at the calculations in a program. If the system does not
include the date or the time you get the value -HUGE(0), that is the
smallest integer number in the model, as output. The vector will
include the following elements: year, month, day, time difference in
minutes, hours, minutes, seconds and milliseconds.
SYSTEM_CLOCK(COUNT, COUNT_RATE, COUNT_MAX)
Subroutine which returns the system time. At least one argument has
to be given. COUNT is a scalar integer which is increased by one for
each cycle up to COUNT_MAX, where it starts once again. If there is no
system clock then -HUGE(0) is returned.
COUNT_RATE is a scalar integer that gives the number of cycles per
second. If there is no system clock the value zero is returned.
COUNT_MAX is a scalar integer which gives the maximum value that
COUNT can reach. If there is no system clock, zero is returned
instead.
----------------------------------------------------------------------
- 71 -
Bit copy routine:
MVBITS(FROM, FROMPOS, LEN, TO, TOPOS)
A subroutine which copies the sequence of bits in position FROMPOS
and has the length LEN to target TO starting in position TOPOS. The
remaining bits are not changed. All quantities have to be integers and
all except TO have to have INTENT(IN) while TO is supposed to have
INTENT(INOUT) and be of the same kind type as FROM. The same variable
can be both FROM and TO. Some natural restrictions apply to the values
of LEN, FROMPOS and TOPOS and you also have to consider the value of
BIT_SIZE.
Random number routines:
A sequence of pseudo random numbers can be generated from a starting
value which is stored as an integer vector. The subroutines offer a
portable interface towards an implementation dependent random number
sequence.
RANDOM_NUMBER(HARVEST)
This subroutine returns in the floating-point number variable HARVEST
one (or several if HARVEST is an array) random numbers between zero and
1.
RANDOM_SEED(size, put, get)
This subroutine resets, or gives information about, the random number
generator. No arguments have to be provided. The output variable SIZE
must be a scalar integer and gives the number of integers (N) the
processor uses for the starting value. The input variable PUT is an
integer vector which puts the starting numbers provided by the user into
the random number generator. The output variable GET (also an integer
vector)reads the present starting value. Example:
CALL RANDOM_SEED ! Initializing
CALL RANDOM SEED (SIZE=K) ! Sets K = N
CALL RANDOM_SEED (PUT = SEED (1:K)) ! Uses the starting value
! given by the user
CALL RANDOM_SEED (GET = OLD(1:K)) ! Returns the present
! starting value
----------------------------------------------------------------------
- 72 -
APPENDIX 6. NAG's Fortran 90
============================
NAGware Fortran 90 compiler
The Numerical Algorithms Group Ltd
Wilkinson House
Jordan Hill Road
Oxford OX2 8DR
infodesk@nag.co.uk Tel: +44 865 511245
Fax: +44 865 310139
NAG offers a Fortran 90 system which behaves as a compiler but in
reality translates to C (except for the PC version). The compiler is
now in version 2.1 for UNIX and 2.01 for MS-DOS. Both work well,
version 2 includes extended error checking both under compilation and
execution. For early experiences please see the article "A First
Encounter With F90" by Michael Metcalf in Fortran Forum, vol. 11, No.
1, pp. 24-32, March 1992. CERN had at that time a total of 80 000 code
lines working with NAG's Fortran 90 compiler, version 1.1. A good
service is obtained from the NAG Response Centre using electronic mail.
The system is available for the following computers: Apollo Domain, DEC
Station Ultrix, DEC VAX/VMS, Hewlett Packard 9000, IBM 6000 and IBM PC,
NeXT, Silicon Graphics, Sun SPARC and Sun 3.
Compilation and execution with the NAG compiler:
You run Fortran 90 in about the same way as you run Fortran 77. The
usual commands are, if the source code "code1.f90" is in the new free
form and source code "code2.f" is in the old fixed form.
% f90 code1.f90 code2.f
and if you have a complete program it can be given an execution name.
% f90 -o program program.f90
Using MS-DOS the compilation command is FTN90 and the source code is
*.F90 for new free form and *.FOR for old fixed form. In order to get
immediate execution you add /LGO as below.
FTN90 PROGRAM.F90 /LGO
When you have the program in more than one file you have to first
compile the various parts separately, and then link with the LINK77
command and obtain an executable program. We now show how to deal with
the two program files PART1.F90 and PART2.FOR and creating a ready to
run program WHOLE.EXE.
FTN90 PART1.F90
FTN90 PART2.FOR
LINK77
$LOAD PART1
$LOAD PART2
$FILE WHOLE.EXE
The symbol $ above is created automatically by the system.
For additional information about f90 under UNIX we recommend the
manual command
----------------------------------------------------------------------
- 73 -
% man f90
I here list the most important switches to the compilation command
-c compilation only (no linking)
-C index control
-fixed the old form of the source code
-free the new form of the source code
-f77 calling convention consistent with the UNIX f77
-g generates information for debugging
-Ipathname Fortran 90 looks automatically for modules in
the present directory, in directories in the
I-list and in /usr/local/lib/f90
-l linking of a library
-Ldir Add dir to the list of directories for library
routines
-o naming the output file
-pg generating execution profile
-S produces source code in the language C
-v give comment about how the compilation is
proceeding
-version gives the compiler's version number
-w suppresses warnings
Peculiarities in the NAG compiler
----------------------------------
o The compiler signals all (?) variables that have been specified but
not used. The compiler therefore performs some data flow analysis.
o The compiler interrupts the compilation if it finds that a
variable that has not been assigned a value is being used.
The compiler thus performs the some data flow analysis.
o The compiler permits Swedish 8-bit characters in character
strings (within apostrophes ' or citation marks ") and in comments.
o The compiler doesn't have support for precision that is larger
than the normal double precision, which means that
SELECTED_REAL_KIND(15) is the highest possible.
o Compiler error messages sometimes state only that an error occurs
on a certain line.
However, it seems that the reported line number is the correct line
number.
o Errors in FORMAT are usually not found at compilation, but only at
execution.
o The maximum unit number UNIT for input/output is 63 instead of
the normal 99. This restriction is removed in version 2.1.
o The last line should be an END with a carriage return (CR),
if not you can get a compilation error! (This has been fixed in
release 2.0, but it is a rather common problem in UNIX).
o With some editors the Tab is stored as a special character, with
some it is automatically replaced with a number of blanks.
Prior to release 2.0 the NAG system did not accept the Tab as
a special character, but it is now an extension.
o Using MS-DOS it is necessary to first start an auxiliary system
with the command DBOS. If you wish to use Windows you have to first
close DBOS with the command KILL_DBOS. If you first start
Windows it is possible to go to the DOS-prompter and start DBOS
there.
---------------------------------------------------------------------
- 74 -
Also at execution of an already compiled program the DBOS
system must be active. It includes among other things a
floating point emulator (used if you do not have floating point
in the hardware).
o The intrinsic function SPREAD can generate too many temporary
arrays, which are not deallocated soon enough (causes memory
shortage).
o The compiler warns for the use of obsolescent concepts, see Appendix 4,
and also for some extensions to Fortran 90.
o The order of the compilation switches is important.
When compiling a program "program.f90" and using modules
in the directory /example/modules and the library "libnag.a"
in the directory /examples/libraries we have to use the
following command
f90 -I/example/modules program.f90 -lnag -L/examples/libraries
System parameters
-----------------
Running the NAG compiler on DEC station ULTRIX, Sun SPARC and IBM PC
(MS-DOS) gives in all cases the following parameters for logical
variables, integers, floating-point numbers and complex numbers. The
identical values are obtained on all the three systems since all these
architectures are based on the IEEE-754. Compare with the numerical
inquiry functions on page 60. PLEASE NOTE THAT NAG HAS SOMEWHAT CHANGED
THE SYSTEM PARAMETERS IN VERSION 2.1. The following values are for the
previous versions.
LOGICAL Default byte word
KIND number = 2 1 2(before rel 2.1)
KIND number = 3 1 3(from rel 2.1)
INTEGER Default int8 int16 int32
KIND number = 3 1 2 3
digits = 31 7 15 31
radix = 2 2 2 2
range = 9 2 4 9
huge = 2147483647 127 32767 2147483647
bit_size = 32 8 16 32
REAL Default single double
KIND number = 1 1 2
digits = 24 24 53
maxexponent = 128 128 1024
minexponent = -125 -125 -1021
precision = 6 6 15
radix = 2 2 2
range = 37 37 307
epsilon = 0.11920929e-06 0.11920929e-06 0.2224460e-15
tiny = 0.11754944E-37 0.11754944E-37 0.22250739-307
huge = 0.34028235E+39 0.34028235E+39 0.17976931+309
COMPLEX Default single double
KIND number = 1 1 2
precision = 6 6 15
range = 37 37 307
---------------------------------------------------------------------
- 75 -
APPENDIX 7. Other Fortran 90 systems
====================================
In 1993 two other Fortran 90 systems were available on the market.
Thinking Machines and Hewlett Packard were very close to full Fortran
90. An interesting product is Forge 90 from Applied Parallel Research,
which translates from Fortran 77 to Fortran 90 and is available for more
UNIX systems.
An interesting study of the three compilers available in 1993 has
been done by John K Prentice, "A Performance Benchmark Study of Fortran
90 Compilers", Fortran Journal, Vol. 5, No. 3, pp. 2 - 7, May/June
1993.
Full Fortran 90 compilera are now available from Cray and Digital
Equipment.
PF90: Portable Fortran 90
-------------------------
Parasoft Corporation
2031 South Myrtle Avenue
Monrovia, CA 91016
f90-info@parasoft.com Tel: +1 818 305 0041
Fax: +1 818 305 9048
Parasoft offers a translation from Fortran 90 (and Connection Machine
Fortran, CM90) to Fortran 77, a library of Fortran 90 internal functions
and a driver to ease the use of the system. Version 2.0 has complete
Fortran 90 and is available for users of Sun SPARC, HP9000, IBM 6000,
IBM 3090 AIX, IBM 9000 ESA and Silicon Graphics.
VAST-90
-------
Pacific-Sierra Research Corporation
Computer Products Group
2901 28th Street
Santa Monica, CA 90405
info@psrv.com Tel: +1 310 314 2338
Fax: +1 310 314 2323
This system offers translation in both directions, from Fortran 90 to
Fortran 77 and from Fortran 77 to Fortran 90. With the translation from
Fortran 77 old constructs are removed, the system also reduces GOTO and
statement numbers, introduces array syntax instead of explicit loops, it
replaces COMMON by module, computed GOTO by CASE and checks possible
execution errors, and a lot more. The new internal functions are,
naturally, included, but some of them are as external functions. The
system is available for Sun SPARC, HP 9000, IBM 6000, Silicon Graphics,
DEC Ultrix, VAX/VMS, Convex and Cray.
The paper by John K Prentice, "Translation of Fortran 77 to Fortran
90 Using VAST-90", Fortran Journal, Vol. 5, No. 4, July/August 1993,
gives a very favourable view of VAST-90 and its performance at
translation. It also includes two actual examples.
---------------------------------------------------------------------
- 76 -
APPENDIX 8. Fortran 90 properties already in Cray CF77
======================================================
Cray has a full Fortran 90 system available, but many of the new
commands in Fortran 90 were already available in Cray Fortran 77, since
Cray Research has been active in the development of the new standard.
The array concepts are almost entirely included, since this area is
close to Cray's general interest in vectorization and parallelization.
See articles in CRAY CHANNELS, Fall 1992, (Vol. 14, No. 4), pp.
18-21.
For the present CF77 version 6.0 the following is already included:
Basic properties:
-----------------
The special characters ", ! and _ are included.
The special character @ is also included but it is used by the system
for special purposes.
The accepted line length can be switched from 72 to 80
characters.
Comments can start with a ! anywhere on the line except
in position 6.
A statement can have up to 99 continuation lines.
A variable name (identifier) can have up to 31 characters, including
the special character _.
Character strings can be delimited with either ' or ".
Octal and hexadecimal constants are available, and
also bit operations according to Appendix 3, item 14.
A repetition factor for new lines has been added, that is
for / in FORMAT.
Only the old form (fix form) of the source code is available.
Control statements:
-------------------
DO, DO WHILE, END DO and WHERE, are available, but not EXIT and
CYCLE.
Specifications:
---------------
IMPLICIT NONE and INCLUDE are available.
The NAMELIST functionality is available, but it can have another
form, in addition to the one from Fortran 90.
The new type of specifications with :: is not available.
File operations.
---------------
All the new specifications to OPEN and INQUIRE are available.
Intrinsic functions:
--------------------
The bit manipulation functions (Appendix 5, item 10) and the bit copy
subroutine MVBITS are available, together with the following array
functions: MERGE, PACK, SPREAD, UNPACK, RESHAPE, SUM, PRODUCT,
MAXVAL, MINVAL, COUNT, ANY, ALL, TRANSPOSE, EOSHIFT, CSHIFT,
MAXLOC, MINLOC, MATMUL and DOT_PRODUCT.
In addition the three bit manipulation functions from HPF are available.
---------------------------------------------------------------------
- 77 -
Arrays:
-------
Adjustable arrays, assumed-size arrays, automatic arrays, array
sections and array assignments are available.
The command ALLOCATABLE is not available.
Additional concepts:
--------------------
Recursion is available with certain special considerations, no RESULT
clause is required.
Pointers are available but they are defined in a different way,
compared with Fortran 90.
You have to be careful with both recursion and pointers on the Cray.
Pointers in Fortran 90 and Cray Fortran are discussed in the article
"Fortran 90 Pointers vs. "Cray" Pointers" by Jeanne Martin in the
Fortran Forum, Vol. 11, No. 2, pp. 17 - 23, June 1992.
Removed extensions:
-------------------
Of compatibility reasons the abbreviation DOUBLE for DOUBLE PRECISION
and the specification DOUBLE COMPLEX were removed. Some system
commands have also been changed.
System parameters on the Cray.
------------------------------
System parameters on the Cray, using the Cray Fortran 90 compiler
are given here. Compare with the NAG table on page 73. Note that Cray
uses the number of bytes for the KIND value. Cray also permits the
values 1, 2 and 4 for logical variables (with the same result as the
default value 8), the values 1 and 2 for integers (with the same result
as the value 4, that is not the same as for the default value 6), and
the value 4 for real and complex variables, (with the same result as
for the default value 8).
LOGICAL Default
KIND number = 8
INTEGER Default int31 int46 int64
KIND number = 6 4 6 8
digits = 46 31 46 63
radix = 2 2 2 2
range = 13 9 13 18
huge = 70368744177663 2147483647 70368744177663 9223372036854775807
bit_size = 46 32 46 46
REAL Default single double
KIND number = 8 8 16
digits = 48 48 96
maxexponent = 8190 8190 8190
minexponent = -8189 -8189 -8189
precision = 14 14 28
radix = 2 2 2
range = 2465 2465 2465
epsilon = 0.71054274E-14 0.71054274E-14 0.25243549E-28
tiny = 0.36672077-2465 0.36672077-2465 0.36672077-2465
huge = 0.27268703+2466 0.27268703+2466 0.27268703+2466
---------------------------------------------------------------------
- 78 -
COMPLEX Default single double
KIND number = 8 8 16
precision = 14 14 28
range = 2465 2465 2465
APPENDIX 9. The historical development of Fortran
==================================================
The following simple program, which uses many different and usual
concepts in programming, is based on "The early development of
Programming Languages" by Donald E. Knuth and Luis Trabb Pardo,
published in "A History of Computing in the Twentieth Century" edited by
N. Metropolis, J. Howlett and Gian-Carlo Cota, Academic Press, New
York, 1980, pp. 197-273. They gave an example in Algol 60 and
translated into some very old languages such as Zuse's Plankalkul,
Goldstine's Flow diagrams, Mauchly's Short Code, Burks' Intermediate PL,
Rutishauser's Klammerausdrucke, Bohm's Formules, Hopper's A-2, Laning
and Zierler's Algebraic interpreter, Backus' FORTRAN 0 and Brooker's
AUTOCODE.
Klammerausdrucke is a German expression, we keep the German
expression also in the Russian and English versions. A direct English
translation is "bracket expression". FORTRAN 0 was really not called
FORTRAN 0, it is just the very first version of Fortran.
The program is given here in Pascal, C and five variants of Fortran.
The purpose of this is to show how Fortran has been developed from a
cryptical, almost machine-dependent language, into a modern structured
high-level programming language.
Program TPK in Pascal for UNIX
------------------------------
program tpk(input,output);
(* Pascal program for Unix *)
var i : integer;
y : real;
a : array [0..10] of real;
function f (t : real) : real;
begin
f := sqrt(abs(t)) + 5*t*t*t
end;
begin
for i := 0 to 10 do read(a[i]);
for i := 10 downto 0 do
begin
y := f(a[i]);
if y > 400 then
writeln(i,' TOO LARGE')
else
writeln(i,y);
end
end.
This program contains variables of the data types integer and
floating-point numbers, and a vector of floating-point numbers. It also
contains internal mathematical functions and a function written by the
user, both forward and backward loops, a conditional statement and
output of both floating-point numbers and characters.
Please note that exponentiation is not included in the Pascal
language, therefore t^3 has to be written t*t*t.
---------------------------------------------------------------------
- 79 -
Program TPK in ANSI C
---------------------
# include
# include
/* Program TPK in ANSI C */
double f (double t);
main()
{
int i;
double y;
double a[11];
for ( i = 0; i <= 10; ++i)
scanf("%lf", &a[i]);
for ( i = 10; i >= 0; i = i -1 )
{
y = f(a[i]);
if ( y > 400 )
{printf("%d",i);
printf(" TOO LARGE\n");}
else
{printf("%d",i);
printf("%lf",y);
printf("\n");}
}
return 0;
}
/* Function */
double f (double t)
{
double temp;
temp = sqrt(fabs(t)) + 5*pow(t,3);
return temp;
}
The language C has the peculiar property that function calls are
normally done in double precision, I have therefore written the whole
program in double precision. This is not really necessary, because
nowadays you can also use the single precision in most realizations.
The exponentiation can be done here with the internal function pow.
[Would you please explain two previous sentences, about single/double?
We are not familiar with the property described.]
[ In the old version of C, calculations were done in double precision,
but you could ask for single, but that was not available when
transferring data between functions and subroutines. See also the
discussion on page 27 of the new revised text. / Bo ]
[Officially the first Fortran version to be spelled Fortran is Fortran 90,
all previous versions should be FORTRAN with upper case. I do not follow
that style, except just a little below.]
---------------------------------------------------------------------
- 80 -
FORTRAN 0
---------
DIMENSION A(11)
READ A
2 DO 3,8,11 J=1,11
3 I=11-J
Y=SQRT(ABS(A(I+1)))+5*A(I+1)**3
IF (400>=Y) 8,4
4 PRINT I,999.
GOTO 2
8 PRINT I,Y
11 STOP
Please note the elegant treatment of the input of the vector A and
that the symbol > (greater than) was available at the beginning of the
Fortran development. Output of text was not available, therefore 999
was used to mark too big numbers. The DO-loop was less elegant, the
digits give starting line and the final line for the loop and where the
execution should be transferred when the loop is terminated. The
conditional statement is also not very user-friendly. Exponentiation is
being done with **.
FORTRAN I
---------
C THE TPK ALGORITHM
C FORTRAN I STYLE
FUNF(T)=SQRTF(ABSF(T))+5.0*T**3
DIMENSION A(11)
1 FORMAT(6F12.4)
READ 1,A
DO 10 J=1,11
I=11-J
Y=FUNF(A(I+1))
IF(400.0-Y)4,8,8
4 PRINT 5,I
5 FORMAT(I10,10H TOO LARGE)
GOTO 10
8 PRINT 9,I,Y
9 FORMAT(I10,F12.7)
10 CONTINUE
STOP 52525
This was the first programming language with comments. The statement
function is used and the backward loop has to be simulated, since a
backward loop (negative index) was not permitted, and neither was index
zero. The conditional statement is the arithmetical one, with jumps to
three different positions depending on whether the arithmetic expression
is less than zero, equal to zero or greater than zero. All function
names have to end with an F. It is possible to print character strings
if you can give the number of characters which will be output (10H in
the case when ten characters are to be outputted). Structured layout
had not yet been invented.
---------------------------------------------------------------------
- 81 -
FORTRAN IV or Fortran 66
-------------------------
C THE TPK ALGORITHM
C FORTRAN IV STYLE
DIMENSION A(11)
FUN(T) = SQRT(ABS(T)) + 5.)*T**3
READ (5,1) A
1 FORMAT(5F10.2)
DO 10 J = 1, 11
I = 11 - J
Y = FUN(A(I+1))
IF (400.0-Y) 4, 8, 8
4 WRITE (6,5) I
5 FORMAT(I10, 10H TOO LARGE)
GO TO 10
8 WRITE(6,9) I, Y
FORMAT(I10, F12.6)
10 CONTINUE
STOP
END
Also here a statement function is being used and the backward loop is
still simulated. The structured layout has been used. The notation
Fortran 66 started to be used in 1978, previously it was just called
Standard FORTRAN or the IBM-notation FORTRAN IV was used also by other
vendors.
Fortran 77
----------
PROGRAM TPK
C THE TPK ALGORITHM
C FORTRAN 77 STYLE
REAL A(0:10)
READ (5,*) A
DO 10 I = 10, 0, -1
Y = FUN(A(I))
IF ( Y . LT. 400) THEN
WRITE(6,9) I,Y
9 FORMAT(I10. F12.6)
ELSE
WRITE (6,5) I
5 FORMAT(I10,' TOO LARGE')
ENDIF
10 CONTINUE
END
---------------------------------------------------------------------
- 82 -
REAL FUNCTION FUN(T)
REAL T
FUN = SQRT(ABS(T)) + 5.0*T**3
END
Also here a statement function ought be used, but is was not popular
at that time, and therefore an external function was being used instead.
The backward loop no longer has to be simulated and the conditional
statement can be given in a more natural way. Structured layout is now
the normal case. The character string can be given within apostrophes
so that you no longer have to count manually the number of characters in
the character string. The statement END in the main program is also
being used for the execution so that STOP is no longer required. Also
the explicit RETURN in the function FUN is now optional.
Fortran 90
----------
PROGRAM TPK
! The TPK Algorithm
! Fortran 90 style
REAL :: Y
REAL, DIMENSION(0:10) :: A
READ (*,*) A
DO I = 10, 0, -1 ! Backwards
Y = FUN(A(I))
IF ( Y < 400.0 ) THEN
WRITE(*,*) I, Y
ELSE
WRITE(*,*) I, ' Too large'
END IF
END DO
CONTAINS ! Local function
FUNCTION FUN(T)
REAL :: FUN
REAL, INTENT(IN) :: T
FUN = SQRT(ABS(T)) + 5.0*T**3
END FUNCTION FUN
END PROGRAM TPK
[Bo, don't you think it would be more feasible to retain unified
I/O formats in all the examples above?]
[Yes, but I am trying to write in the style of the period,
not in the best possible way. / Bo]
Instead of an external function a local function is now used. The
backward loop has not to be simulated. The conditional statement can be
given in a more natural way and the symbol < (less than) has returned
from Fortran 0. The specifications have a new appearance. Structured
layout is used always and all statement numbers have been removed. The
character string is given within apostrophes and the standard formats
and standard units are are used. Comments can be given on the line.
-----------------------------------------------------------------------------
- 83 -
APPENDIX 10. High Performance Fortran
=====================================
Very important, at the introduction of a new programming language
standard, is nowadays that the language should permit efficient
compilation and execution not only on conventional computers, but also
on parallel and massively parallel systems. Somewhat simplified it
could be stated that Fortran 90 is efficient on conventional computers
and on vector processors, but less efficient on parallel processors.
The reason for this is that in Fortran it does not exist any commands
for partition into parallel processes, except for those that follow
directly from the array concept, with its operators and functions.
A group based in Houston, High Performance Fortran Forum, has
developed a proposal in the form of an extension to Fortran 90. The aim
of this project High Performance Fortran or HPF is to offer a portable
language which permits an efficient use on different parallel systems.
The project has issued a final proposal in May 1993, and aims towards a
de facto standard, and not formal standard from ANSI or ISO. In order
to facilitate the introduction and general acceptance of HPF the group
has also defined a subset based on Fortran 77 and just a few parts of
Fortran 90. A number of manufacturers were involved in the group, and
the proposal will hopefully get a fast acceptance, with many
implementations. Those parts of the proposal that were controversial,
and therefore requires further study, is available in a separate
document "Journal of Development".
The proposal is a Fortran-based language with facilities to control
distribution of arrays onto distributed memory parallel computers and
includes extensions for
Parallel execution
Good performance both on MIMD and SIMD computers
Distribution of data on the available processors
Location of data within a specific processor
The complete HPF is based on Fortran 90 and includes
New directives
New syntax (the statement FORALL and new intrinsic functions)
Restrictions in the rules for storage
The Subset HPF is instead based on Fortran 77 in order to become
implemented already in 1993 and includes
MIL-STD-1753 (US military extension to Fortran 77, it is also
included in Fortran 90)
Array-operations
All HPF-constructs, but with certain limitations in their
functionality
The directives have been given such a form that it will be possible
to use them as ordinary statements in a future version of Fortran. At
present they will have to be written as comments.
The first HPF system seems to be from Digital Equipment.
CHPF$ directive ! Fortran 77 and fix form Fortran 90
!HPF$ directive ! Fortran 90 (both fix and free form)
-----------------------------------------------------------------------
- 84 -
Data storage
------------
HPF has directives on how data ought to be stored. These directives
do not influence the computational results, but they may influence the
computational speed (hopefully increasing it). The reason for
introducing the directives for data storage are the following two very
simple observations:
o An operation on several data objects ought to be faster if all
the objects are on the same processor.
o It ought to be possible to perform many operations at the same
time if the operations are performed on different processors.
Using the data storage directives we can distribute the data on the
various processors. Sometimes it is desirable first to align an array
against a template, and then distribute the template on the processors.
A first example can explain the situation.
REAL A(1000), B(1000), C(1000), X(500), Y(0:501)
INTEGER INX(1000)
!HPF$ PROCESSORS PROCS(10)
!HPF$ DISTRIBUTE A(BLOCK), B(BLOCK) ON PROCS
!HPF$ DISTRIBUTE INX(BLOCK) ON PROCS
!HPF$ DISTRIBUTE C(CYCLIC) ON PROCS
!HPF$ ALIGN X(I) WITH Y(I-1)
A(I) = B(I) ! (1)
X(I) = Y(I-1) ! (2)
A(I) = C(I) ! (3)
A(I) = A(I-1) + A(I) + A(I+1) ! (4)
C(I) = C(I-1) + C(I) + C(I+1) ! (5)
X(I) = Y(I) ! (6)
A(I) = A(INX(I)) + B(INX(I)) ! (7)
We here work with 10 processors, and we have distributed the floating
point vectors A and B and also the integer vector INX block wise over
these processors. We then have the first one hundred elements from each
one of the vectors on the first processor, the next one hundred on the
next processor, and so on. The floating point vector C is however
distributed cyclically, the first, eleventh, twenty-first etc. elements
are on the first processor, the second, twelfth, twenty-second etc.
elements are on the second, and so on.
Without instructing the system exactly how the two vectors X and Y
are to be distributed, we require that the I-th element of the vector X
is on the same processor as the (I-1)-th element of the vector Y.
The assignments in the cases (1) and (2) above can thus be done
without any inter-processor communication, since all data all the time
is on the same processor.
The assignment in case (3) will in most instances be between
different processors, only in exceptional cases will the data be on the
same processor.
The assignment in case (4) on the other hand will in most instances
be between data on the same processor, except in the exceptional cases
at block boundaries.
-----------------------------------------------------------------------
- 85 -
The case (5) looks like case (4), but in this case the data is all
the time on different processors.
For case (6) we only know that X(I) and Y(I-1) are on the same
processor, but this gives no information whatsoever on how X(I) and Y(I)
are distributed. At block wise storage both of them are normally on the
same processor, at cyclic storage they are always on different processor
(provided that the cycles are the same).
We do not know very much about the distribution in case (7). It is
however easy to realize that A(INX(I)) and B(INX(I)) always are on the
same processor, and that also A(I) and INX(I) are on the same processor.
It is however not likely that these two sets will coincide.
Also the formal arguments in a function or a subroutine can be stored
in different ways using directives. Assume that we have a subroutine
with two arguments, both floating point vectors.
SUBROUTINE CAROLUS (KARL, CHARLES)
REAL, DIMENSION (1718) :: KARL, CHARLES
!HPF$ ALIGN WITH * :: KARL
!HPF$ ALIGN WITH KARL :: CHARLES
...
END
This is called with
REAL, DIMENSION (1718) :: GUSTAV, ADOLF
...
CALL CAROLUS (GUSTAV, ADOLF)
This means that the first formal argument KARL is distributed in the
same way as the actual argument GUSTAV. The second formal argument
CHARLES is distributed in the same way as the first formal argument, in
this case in the same way as GUSTAV, and not as the second actual
argument ADOLF.
Another example is when you require the four corner submatrices of
the size N*N from a larger matrix of the size (N+1)*(N+1). In this case
we call the large matrix EARTH and the submatrices are called NW, NE,
SW, and SE. In some cases the large matrix is not required to be used
in any calculations, its only purpose is to assist at the distribution
of the data in the submatrices. It is therefore a TEMPLATE, which is
not stored.
!HPF$ TEMPLATE, DISTRIBUTE (BLOCK, BLOCK) :: EARTH (N+1, N+1)
REAL, DIMENSION (N, N) :: NW, NE, SW, SE
!HPF$ ALIGN NW(I, J) WITH EARTH(I, J)
!HPF$ ALIGN NE(I, J) WITH EARTH(I, J+1)
!HPF$ ALIGN SW(I, J) WITH EARTH(I+1, J)
!HPF$ ALIGN SE(I, J) WITH EARTH(I+1, J+1)
Since a TEMPLATE does not represent any real storage, it can not be
part of any COMMON. The alignment subscripts can be just a little more
complicated than the old array indices in Fortran 66, they are of the
general form m*i+n, where i is a variable which is permitted to appear
only once in the expression, and m and n are constants (PARAMETER or
numerical constant).
A somewhat more complicated directive is
!HPF$ ALIGN A(:) WITH D(:,*)
-----------------------------------------------------------------------
- 84 -
which means that a copy of the vector A is associated with each
column of the matrix C. An example on this follows, where it is also
shown that the HPF syntax permits several ways for expressing identical
requirements.
!HPF$ TEMPLATE, D1(N), D2(N, N)
REAL, DIMENSION (N, N) :: X, A, B, C, AR1, AR2, P, Q, R, S
!HPF$ ALIGN X(:,*) WITH D1(:)
!HPF$ ALIGN (:,*) WITH D1 :: A, B, C, AR1, AR2
!HPF$ ALIGN WITH D2, DYNAMIC :: P, Q, R, S
The following is a more complete example, where an intrinsic function
finds the present number of processors, which is used at the
distribution of the arrays. In addition a usual (external) function is
being used.
PARAMETER ( N = NUMBER_OF_PROCESSORS() )
!HPF$ PROCESSORS MPP(N)
!HPF$ TEMPLATE T(1000), S(5000)
!HPF$ DISTRIBUTE T(BLOCK) ONTO MPP ! Block size may
!HPF$ DISTRIBUTE S(CYCLIC) ONTO MPP ! be specified
REAL, DIMENSION (1000) :: X, Y, Z
REAL, DIMENSION (5000) :: V, W
!HPF$ ALIGN WITH T :: X, Y, Z
!HPF$ ALIGN WITH S :: V, W
REAL, DIMENSION (10,1000) :: A
!HPF$ ALIGN WITH T :: A(*,:) ! Vector of columns
...
TEMP = F(V)
...
END
FUNCTION F(X)
REAL, DIMENSION (:) :: X
!HPF$ ALIGN WITH * :: X ! Distribute the formal argument
! X as the real argument V
REAL, DIMENSION (SIZE(X)) :: S
!HPF$ ALIGN WITH X :: S ! Distribute the local vector S as
! the formal argument X, or as the
! actual argument V
...
RETURN
END
You can specify different distributions along different dimensions.
The following specifications
-----------------------------------------------------------------------
- 87 -
REAL A(100,100), B(100,100), C(200)
!HPF$ DISTRIBUTE A(BLOCK,*), B(*,CYCLIC), C(BLOCK(5))
means that the first processor of a four processor computer stores the
following array sections
A(1:25, 1:100)
B(1:100, 1:97:4)
C(1:5), C(21:25), C(41:45), C(61,65), C(81:85),
C(101:105), C(121:125), C(141:145), C(161,165),
C(181:185)
It is also possible to distribute several dimensions in completely
independent ways,
REAL D(8,100,100)
!HPF$ DISTRIBUTE D(*,BLOCK, CYCLIC)
means that the first processor of a four processor computer,
configured as a 2*2 matrix, stores the following array sections
D(1:8, 1:50, 1:99:2)
In addition to the static directives discussed above, there are also
the two dynamic directives REDISTRIBUTE and REALIGN, who permit an array
to switch its distribution within the subroutine. At the use of a
subroutine it does exist three different possibilities for the formal
arguments.
o A formal argument can inherit the distribution from the actual
argument.
o A formal argument can be redistributed to a requested
distribution.
o A formal argument can be specified to have a certain distribution,
and is an error if the actual argument does not have that
distribution.
At the exit from a subroutine it is necessary that the effect of any
redistribution is removed.
Execution
---------
A new statement FORALL is introduced as an alternative to the
DO-loop. The intent with the FORALL statement is that its content can
be executed in any order, independent of the index! It therefore gives
the possibility of a parallel implementation. Some rules are however
introduced in order to avoid side effects. Some simple examples on
FORALL follow.
FORALL (I = 1:N, J = 1:N) H(I, J) = 1.0/REAL(I+J-1)
FORALL (I = 1:N, J = 1:N, Y(I, J) .NE. 0.0) &
X(I, J) = 1.0/Y(I, J)
FORALL (I = 1:N) A(I, I+1:N) = 3.141592654
The first of these define a Hilbert matrix of order N, the second
inverts the elements of a matrix, avoiding division with zero. In the
third example all elements above the main diagonal of the matrix are
assigned the value of the mathematical constant pi.
-----------------------------------------------------------------------
- 88 -
In all the three statements above, FORALL can be considered as a
double loop, which can be executed in arbitrary order. The general form
of the FORALL statement is
FORALL ( v1 = l1:u1:s1, ... , vn = ln:un:sn, mask ) &
a(e1, ... , em) = right_hand_side
and is evaluated according to certain well specified rules, in
principle all indices are evaluated first.
In addition there is a FORALL construct, which has to be
distinguished from the single line FORALL statement. The construct is
interpreted in the same way as if its internal statements were written
as FORALL statements in the same order. This restriction is important
in order to achieve a unique computational result. Within a FORALL
construct calls of functions or subroutines are permitted only if these
subprograms are pure (without side effects). First we give a simple
example of a FORALL construct.
REAL, DIMENSION(N, N) :: A, B
...
FORALL (I = 2:N-1, J = 2:N-1)
A(I, J) = 0.25*(A(I, J-1)+A(I, J+1)+A(I-1, J)+A(I+1, J))
B(I, J) = A(I, J)
END FORALL
When these statements have been executed the arrays A and B have
identical values in the internal points, while B has kept its previous
values on the boundaries.
In addition a directive INDEPENDENT has been introduced for both DO
loops and FORALL constructs. This directive is placed immediately
before the DO statement or FORALL construct, and is valid until the
corresponding END DO (or the old form of terminating a DO loop) or END
FORALL. The directive assures the system that this part of the program
can be executed in an arbitrary order, including parallel, without any
computational differences in the result (no semantic change).
In the example below it is thus assured that the integer vector P
does not have any repeated values (which would have meant that last one
wins at a normal sequential execution). A potential conflict at
parallel execution is thus avoided. It is also implicitly assured that
all values of P are within the permitted bounds of 1 and 200.
REAL, DIMENSION(200) :: A
REAL, DIMENSION(100) :: B
INTEGER, DIMENSION(100) :: P
...
!HPF$ INDEPENDENT
DO I = 1, 100
A(P(I)) = B(I)
END DO
It is also possible to indicate that certain parts of a nested loop
shall be considered as independent. In the example below the innermost
loop is not independent since each element of A is assigned several
times, for all values of I4.
-----------------------------------------------------------------------
- 89 -
REAL, DIMENSION(N, N, N) :: A, B, C
...
!HPF$ INDEPENDENT, NEW (I2)
DO I1 = 1, N1
!HPF$ INDEPENDENT, NEW (I3)
DO I2 = 1, N2
!HPF$ INDEPENDENT, NEW (I4)
DO I3 = 1, N3
DO I4 = 1, N4 ! The innermost loop is NOT
! independent !
A(I1, I2, I3) = A(I1, I2, I3) &
+ B(I1, I2, I4) * C(I2, I3, I4)
END DO
END DO
END DO
END DO
The NEW clauses are required, since the inner loop indices are
assigned and used in different iterations of the outer loops.
A DO loop and a FORALL construct with the same content are equivalent
if they both have been given the INDEPENDENT directive.
Intrinsic functions
-------------------
In addition to all the intrinsic functions of Fortran 90 the language
extension HPF introduces a few new, and modifies two of those already in
Fortran 90.
To the functions MAXLOC and MINLOC have been added an optional
parameter DIM, which informs along which dimension the maximum or
minimum is to be taken. An example illustrates how it works on MAXLOC.
Without the extra parameter the index for the largest element is
returned. If the optional parameter DIM = 1 the row numbers for the
largest element in each column are returned. If the optional parameter
DIM = 2 the column numbers for the largest element in each row are
returned. The corresponding applies to MINLOC. These new facilities
are the same as those that already exist in Fortran 90 for MAXVAL and
MINVAL.
0 -5 8 -3
A = 3 4 -1 2
1 5 6 -4
gives the following values
MAXLOC(A) = (/ 1, 3 /)
MAXLOC(A, DIM = 1) = (/ 2, 3, 1, 2 /)
MAXLOC(A, DIM = 2) = (/ 3, 2, 3 /)
The following completely new functions have been added. The inquiry
functions are to be intrinsic, but the others may instead be available
in a library as external functions.
-----------------------------------------------------------------------
- 90 -
Inquiry functions:
NUMBER_OF_PROCESSORS(DIM) Total number of processors
PROCESSORS_SHAPE>() Shape for the processors
Example. For a DECmpp 12000 Model 8B we get
NUMBER_OF_PROCESSORS() 8192
NUMBER_OF_PROCESSORS(DIM=1) 128
NUMBER_OF_PROCESSORS(DIM=2) 64
PROCESSORS_SHAPE() (/ 128, 64 /)
while on an ordinary workstation we get
NUMBER_OF_PROCESSORS() 1
PROCESSORS_SHAPE() (/ 1 /)
Bit Manipulation Functions: These three use the model in section
13.5.7 of the Fortran 90 standard. Note that the results are machine
dependent, since they require the number of bits in an integer, which
in Fortran 90 is available with the intrinsic function BIT_SIZE.
POPCNT Gives the number of 1-bits in the integer argument
POPPAR Gives the result 1 if the number of 1-bits in the
integer argument is odd, else the result 0.
LEADZ Gives the number of leading zeros in the integer
argument.
Other functions:
The function ILEN is an integer function which returns the number of
bits required to store a signed integer in the form of a 2 complement.
Examples of its use are that 2**ILEN(N-1) rounds an arbitrary positive
integer N upwards towards the closest power of 2, and 2**(ILEN(N)-1)
rounds down.
The following routines are not required to be intrinsic but may be in
a normal external library, and in that case they require a USE
statement. They are not explained below, we refer the reader to the
HPF standard.
The array reduction functions IALL, IANY, IPARITY and PARITY are
available and they correspond to the following intrinsic functions of
Fortran 90, namely IAND, IOR, IEOR and the operator .NEQV.
A large number of functions are available to gather and scatter data,
they have names of the form XXX_SCATTER, where XXX can be SUM, COUNT,
PRODUCT, MAXVAL, MINVAL, IALL, IANY, IPARITY, ALL, ANY and PARITY.
For parallel operations there are the functions XXX_PREFIX and
XXX_SUFFIX, where XXX has the same possibilities as for XXX_SCATTER.
Example: SUM_PREFIX sums successively the elements of the array, the
first remains unchanged, the second becomes the sum of the first two,
and so on. With SUFFIX the summation is done in the other direction
(backwards).
In addition there are two sorting functions, GRADE_UP and GRADE_DOWN.
Operations for parallel input and output are being considered, and can
be found in the Journal of Development.
-----------------------------------------------------------------------
- 91 -
Storage rules
-------------
For Fortran 90 there exist certain rules about storage of variables,
especially for arrays who have to be stored in a well specified order,
with the first index varying fastest. There are also special rules for
the storage of variables that are in COMMON and/or EQUIVALENCE. These
rules on sequential storage can conflict with the distribution requested
by the HPF statements ALIGN, REALIGN, DISTRIBUTE and REDISTRIBUTE. The
Fortran model works with linear storage with one storage unit for
integers and ordinary floating point variables, and two storage units
for double precision floating point numbers.
Fortran has storage association and sequence association. The
Fortran 90 standard states in (14.6.3) that storage association is the
association of two or more data objects that occurs when two or more
storage sequences share or are aligned with one or more storage units,
and in (12.4.1.4) that sequence association is the order that Fortran
requires when an array is associated a formal argument. The rank and
shape of the actual argument neeed not agree with the rank and shape of
the dummy argument, but the number of elements in the dummy argument
must not exceed the number of elements in the element sequence of the
actual argument. If the dummy argument is assumed size, the number of
elements in the dummy argument is exactly the number of elements in the
element sequence.
Note that HPF has no problem with array parameters distributed over
the processors, as long as both the actual and the dummy arguments have
the same rank and shape. It is when the properties of Fortran, with
respect to COMMON and EQUIVALENCE, are used too much, that we get into
problems. If we use a subroutine that contains the following
specifications
SUBROUTINE HOME(X)
DIMENSION X(20,10)
it can be called with CALL HOME(Y(2,1)) provided that the array X is
specified SEQUENTIAL in the subroutine HOME and the array Y is also
specified SEQUENTIAL in the calling program unit. The simple call CALL
HOME(Y) is permitted if X and Y are both sequential, or if X and Y are
dimensioned in exactly the same way.
A directive SEQUENCE has therefore been introduced in order to inform
the system that a variable or a COMMON block shall be considered
sequential. Using HPF the normal case is that COMMON blocks are
considered non-sequential, while variables normally also are
non-sequential (except in a few cases, as when the variable is part of a
sequential COMMON block or it is an array with assumed dimension).
In order to get old programs to produce the same results with HPF the
following is recommended:
o Include a facility to check that all parts of COMMON blocks
agree within the whole program.
o Include a facility so that all COMMON blocks are considered
sequential, if no explicit directive for the opposite case
has been given.
-----------------------------------------------------------------------
- 92 -
HPF Subset
----------
The following news from Fortran 90 are included in the HPF Subset:
Ordinary concepts:
Identifiers with up to 31 characters
including the underline symbol _
Comments starting with an exclamation mark !
Array concepts
Array sections
Array assignment
The statements WHERE and ELSEWHERE
Array valued external functions
Automatic arrays
Allocatable arrays
Arrays with assumed shape
New intrinsics
Array valued new intrinsics
Additional concepts from Fortran 90
Optional arguments
Keyword arguments
INTERFACE (but not generic or for modules)
Type specifications (but not KIND and TYPE)
In addition the military extension MIL-STD-1753 to Fortran 77 is
included. It has DO WHILE, END DO, IMPLICIT NONE and INCLUDE; the bit
manipulation functions BTEST, IAND, IBCLR, IBITS, IBSET, IEOR, IOR,
ISHFT, ISHFTC and NOT; the bit copy routine MVBITS; and the binary,
octal and hexadecimal constants.
The following from the complete HPF is however not included in the
HPF Subset:
The directives DYNAMIC, EXTRINSIC, LOCAL, PROCESSOR_VIEW,
PURE, REALIGN and REDISTRIBUTE
The construct FORALL ... END FORALL
The HPF-library
The HPF_LIBRARY module
In addition there are certain restrictions in the subset on the
permitted values of the argument DIM to MAXLOC and MINLOC, and there is
no definition of extrinsic subprograms and no binding to Fortran 90.
------------------------------------------------------------------------------
- 93 -
APPENDIX 11. Explanation of certain terms
=========================================
[ Listed alphabetically according to the English word, therefore the
actual page number does not agree with the Swedish version]
adjustable array The array in a subprogram has a variable
dimension just
as in Fortran 77 through that both the
array name and the dimension (all dimensions)
are dummy arguments. The actual shape is
therefore passed from the actual arguments.
SUBROUTINE sub (a, na)
REAL, DIMENSION (na) :: a
allocatable array An array specified as ALLOCATABLE with
a certain type and rank. It can later
be allocated a certain extent with
the ALLOCATE statement. The corresponding
storage area can released with the
DEALLOCATE statement.
array Dimensioned quantity (variable or constant).
array assignment Assignment of a whole array to another one
is permitted provided they have the same
shape, and also if the expression on the right
hand side is a scalar, in which case all
the elements on the left hand side will be
assigned this same value.
Array sections can be used here instead
of arrays.
array function A function which operates on an array and
returns an array or a scalar, or a function
which operates on a scalar and returns
an array.
array section Part of an array (the part can have a rank
greater than one, but array sections with
rank one are the most common cases).
assumed-shape array An array in a subprogram which is not local
and has a certain data type and a certain
rank.
SUBROUTINE sub (a)
REAL, DIMENSION (:, :, :) :: a
The extent is determined from an explicit
INTERFACE in the calling program unit,
and an explicit specification of the shape
of A there (including the extent).
INTERFACE
SUBROUTINE SUB (A)
REAL, DIMENSION (:, :, :) :: A
END SUBROUTINE SUB
END INTERFACE
--------------------------------------------------------------------------------
- 94 -
assumed-size array The array has a variable dimension, like in
Fortran 77, through that both the array name
and the parameters of the dimension are
dummy arguments, except the last dimension,
which is given by a *
SUBROUTINE sub (a, na1, na2)
REAL, DIMENSION (na1, na2,*) :: a
attribute At the specification of a variable the
additional properties PARAMETER,
PUBLIC, PRIVATE, INTENT, DIMENSION,
SAVE, OPTIONAL, POINTER, TARGET,
ALLOCATABLE, EXTERNAL and INTRINSIC
are called attributes. They can also
be given by statements.
automatic array Array in a subprogram where the array is
local but the parameters for the dimensions
are among the dummy arguments.
SUBROUTINE sub (i, j, k)
REAL, DIMENSION (i, j, k) :: x
BLOCK DATA A BLOCK DATA program unit contains
definitions (initial values) which are
to be used in the other program units.
BLOCK DATA is now being replaced with the
more powerful module.
extent The number of elements along the
various dimensions.
rank and its extent.
EXTERNAL The functions and subroutines you write
or obtain from application libraries.
The opposite is INTRINSIC.
When external functions or subroutines
are used as actual parameters at the reference
of another subprogram they have to be
specified as EXTERNAL.
function A program unit which returns a function value
(or several, if array valued) in its name
and has a number of arguments. Functions
without arguments are permitted, it is also
permitted but not advisable that a
function changes any of its arguments.
generic A generic function can have the property
that the data type of the argument gives
the same data type to the function. The value
SIN(1.0D0) is also more accurate (calculated
to a higher precision) than that of SIN(1.0).
In addition, there exist intrinsic functions,
for example REAL() which can have arguments
of various data types, but always returns a
value of a certain data type. For generic
subroutines only the arguments are adjusted.
------------------------------------------------------------------------------
- 95 -
interface Since all the different program units in a
Fortran program are compiled completely
independent all information on the arguments
have to be transferred manually. In Fortran 77
this was done with very long and cumbersome
argument lists. In Fortran 90 an INTERFACE
can be used. This has to be used in the
following cases
a) modules that use user defined data types
b) call with keyword arguments or optional arguments
c) user defined generic routines
d) use of assumed-shape arrays
e) use of arrays specified with pointers
f) at the definition of a new interpretation of an OPERATOR
g) to increase the chance that INTENT has any effect
(implementation dependent case)
h) call of an array valued function
i) If IMPLICIT NONE is being used, an INTERFACE may be
required at the use of a subprogram as the actual argument
internal function A function which is local to a certain
program unit, and is after a CONTAINS
statement. All the variables in the
program unit are directly available.
An internal function can not be used
outside the program unit where it is defined.
This can be useful in order to avoid
name conflicts between functions and
subroutines from different libraries.
An internal function can not be used as
an argument! It also exists internal
subroutines.
intrinsic Those functions and subroutines that are
delivered together with the compiler are
called intrinsic, and have some special
properties, they do not usually have to
be specified. Since Fortran 90 has so
many intrinsic functions, some manufacturers
may choose to let some of them instead be
available in a usual library, and thus be
external.
The opposite to intrinsic is external.
If intrinsic functions or subroutines
are used as arguments at the call of
some other subprogram they must be
specified as INTRINSIC. Note the very
important difference between an intrinsic
and an internal subprogram.
------------------------------------------------------------------------------
- 96 -
loop With a loop we mean a number of executable
statements that are repeated a number of
times,
the usual case is a DO-loop, but also an
IF-statement (arithmetic or logical or the
newer IF...THEN...ELSE...ENDIF or a CASE
construct can be called a loop.
main program Each Fortran program has to consist of
exactly one main program and and any number
of subroutines, functions, modules
and BLOCK DATA program units.
A main program may start with the statement
PROGRAM name
and must be terminated with a statement
END PROGRAM name
or the simpler
END
module A module contains specifications and
definitions to be used in other program
units. Replaces BLOCK DATA.
nested Inside a DO-loop can be further DO-loops,
inside an IF...THEN...ELSE...ENDIF or
CASE can also be additional constructs.
Sometimes also IF and CASE constructs
are called loops.
program unit The collective reference to a main program,
a subroutine, a function, a module or a
BLOCK DATA program unit.
rank The number of dimensions.
WARNING: Not the mathematical rank.
recursive A function or a subroutine that calls itself.
Permitted from Fortran 90.
result variable A function is called by its function name,
and on return that name contains the value
of the function. If the function is recursive
a special result variable has to be used
inside the function in order to store the
result temporarily, but it is still returned
in the usual way via the function name.
The result variable is also useful for
array functions, even in the non-recursive case.
It is always permitted in a function
specification to use the result variable.
shape The shape of an array consists of its
rank and extent.
statement function A function local to a program unit, is
located between the ordinary specifications
and the executable statements. Very simple,
available already in FORTRAN I. A more general
concept is the internal function after a
CONTAINS statement. A statement function may
not be used as an argument!
------------------------------------------------------------------------------
- 97 -
subroutine A program unit which does not return any
value through its name and has a number
of arguments. It is permitted for a
subroutine to change the values of its
arguments. A subroutine can however have
other tasks to perform, quite often related
to input and output.
A subroutine does not have a data type.
APPENDIX 12. Argument association
=================================
Basic concepts.
---------------
Dummy argument is a quantity in the specification of a procedure,
function or subroutine.
Actual argument is a quantity in the call of a procedure, function or
subroutine. The call is made from the calling procedure or program
unit, often from the main program.
Pascal.
-------
When an argument is not specied as var then the content of the
quantity will be copied to a new storage area, which will be used within
the procedure. This means that , if the quantity is changed within the
procedure, the new value will not be returned to the calling procedure.
Such an argument is said to be a value parameter. With the terminology
of Fortran 90 it should have the attribute INTENT(IN).
In Pascal the assignment of a new value to a variable not specified
var is completely legal, but the new value is not returned to the
calling program.
When an argument is specied as var then both the dummy and actual
arguments will refer to the same storage area. This means that, if the
quantity is changed within the procedure, the new value will be returned
to the calling procedure. Such an argument is said to be a variable
parameter. With the terminology of Fortran 90 it should have the
attribute INTENT(INOUT).
In the following very small program the procedure sub has x and y as
value parameters, and z as a variable parameter. Here the dummy
arguments are x, y and z, while the actual arguments are a, b and c.
program example;
var a, b, c : real;
procedure sub(x, y : real; var z : real);
begin
y := y*y;
z := x + y;
end;
begin
a := 1;
b := 10;
sub(a, b, c);
writeln(' a = ', a); (* Becomes 1 *)
writeln(' b = ', b); (* Becomes 10 *)
writeln(' c = ', c); (* Becomes 101 *)
end.
-----------------------------------------------------------------------------
- 98 -
Fortran.
--------
A dummy argument, which in Fortran 90 has been specified as
INTENT(IN) in the called program unit, should generate an error message
if the argument is assigned a new value within the called program unit.
In the following small Fortran program the subroutine SUB has X as an
input variable, Y as an input and output variable, and Z as an output
variable. The intent specification should therefore be as below.
Please note that use of INTENT is optional but highly recommended. All
compilers do not yet support it fully (all possible checks are not
performed. Note that we have chosen different names for the actual
arguments and for the formal arguments in the INTERFACE and in the
SUBROUTINE. Usually the same name is chosen for all three cases.
Here the dummy arguments are X, Y and Z, the actual arguments are A,
B and C, and the dummy arguments in the INTERFACE are S, T and U.
PROGRAM EXAMPLE
IMPLICIT NONE
INTERFACE
SUBROUTINE SUB(S, T, U)
REAL, INTENT(IN) :: S
REAL, INTENT(INOUT) :: T
REAL, INTENT(OUT) :: U
END SUBROUTINE SUB
END INTERFACE
REAL :: A, B, C
A = 1
B = 10
CALL SUB(A, B, C)
WRITE(*,*) ' A = ', A ! Becomes 1
WRITE(*,*) ' B = ', B ! Becomes 100
WRITE(*,*) ' C = ', C ! Becomes 101
END PROGRAM EXAMPLE
SUBROUTINE SUB(X, Y, Z)
IMPLICIT NONE
REAL, INTENT(IN) :: X
REAL, INTENT(INOUT) :: Y
REAL, INTENT(OUT) :: Z
Y = Y**2
Z = X + Y
END SUBROUTINE SUB
In Fortran there is really only one kind of argument usage, called
argument association. There may still however be two cases two
consider. If the argument is a scalar then the value of the actual
argument may be copied to a memory location for the dummy argument, and
on return (provided the value has been changed) the copy process will be
repeated in the opposite direction, from the called program unit to the
to the calling
------------------------------------------------------------------------------
- 99 -
program unit or the main program. Please note that the first copy
might be done to a register and not necessarily to the main RAM memory.
The copy process back to the calling program unit naturally fails if the
actual argument is a constant or an expression , for example 249 or
SIN(X) or X+2, instead of a variable name. Normally no error message is
obtained. The intent of the INTENT attribute is to give the compiler
greater possibilities to check against similar inconsistencies.
Also note that if F(X) is a function which is called, then A in F(A)
is a variable name, but (A) in F((A)) is a more general expression.
This possibility to introduce an extra set of parenthesis in order to
move from a variable to an expression works sometimes, but is very
implementation dependent, and therefore strongly discouraged.
If the formal argument is an array then the mentioned copy process
could mean a very large memory requirement, because the arrays would be
stored twice, both in the calling and the called program unit. This is
solved in such a way that the actual argument is used directly in the
called subprogram. The actual array argument must have the same data
type, kind and shape as the dummy argument, but there are a few
exceptions from this general rule. Fortran assumes that the actual
array has all the properties that were specified for the dummy array
argument. The position for each array element is therefore calculated
from this assumption. If the actual array is too different from the
dummy array the wrong element can be referenced, with even fatal
consequences. Usually this rule can be utilized in a satisfactory
manner.
If the arguments are pointers or names of functions or subroutines
the discussion above is not completely valid.
In the good old days IBM referred to that a variable was received by
value regarding what was discussed above for scalars, and received by
location regarding arrays (1977 about FORTRAN IV for IBM 360, an
extension of Fortran 66). There even was a method to force received by
location. The difference was rather important, but that is not the case
any more. We used to have some problems with synchronization between
index variables and variables in COMMON.
Storage rules.
--------------
In Fortran there are certain rules about storage of variables,
especially for arrays which have to be stored in a well defined order
with the first index varying fastest, and also for variables in COMMON
and EQUIVALENCE. The Fortran model is completely linear, with a certain
memory requirement for integers and for floating point values in single
precision, and the double requirement for floating point values in
double precision.
Fortran has storage association and sequence association. The
Fortran 90 standard states in (14.6.3) that storage association is the
association of two or more data objects that occurs when two or more
storage sequences share or are aligned with one or more storage units,
and in (12.4.1.4) that sequence association is the order that Fortran
requires when an array is associated a formal argument. The rank and
shape of the actual argument need not agree with the rank and shape of
the dummy argument.
----------------------------------------------------------------------
- 100 -
INDEX
=====
In this index I have indicated pages with an explanation of a certain
term with a boldface font, and pages where there are only an occurrence
of the term with the ordinary font. This separation is however not
included in this ASCII-version.
The FORMAT-letters are not included below, for those I refer the
reader to Appendix 2, pp. 46-48.
[ THE ENGLISH INDEX IS IN A PRELIMINARY VERSION ]
[ ab Seitenzahl 71 ist 1, ab Seitenzahl 81 ist 2 zu addieren ]
A
ABS 48, 78-80
ACCESS 48
ACHAR 58
ACOS 48
ACTION 48
adjustable array 20, 92
ADJUSTL 58
ADJUSTR 58
ADVANCE 25, 48
AIMAG 48
AINT 48, 58
Algol 49, 76
align 82-83
ALIGN 82-84, 89
ALL 26, 62, 75, 88
Alliant 74
ANINT 48, 58
ANY 62, 75, 88
ALLOCATABLE 53, 92
allocatable array 20, 53, 90-91
ALLOCATE 21, 22, 37, 53
ALLOCATED 63
alternative representation 49
Apollo 71
APOSTROPHE 48
APPEND 48
argument (parameter) 95-97
arithmetic IF-statement 44, 57
array 16, 19-20, 52, 81, 89-92
array section 19, 36, 90-91
ASA 48
ASIN 48
ASIS 48
ASSIGN 44, 57
assigned GOTO-statement 10, 44, 57
ASSOCIATED 21, 69
assumed dimension 20, 75, 89, 91
assumed shape 20, 35, 54, 91
ATAN 7, 29, 48
ATAN2 48
attribute 11, 14, 21, 27, 53-54, 56, 92
automatic array 20, 53, 90-91
-----------------------------------------------------------------------
- 101 -
B
BACKSPACE 46
backward compatibility 2, 57, 89
binary number 26, 47, 55, 60, 90
bit 55, 60, 75, 88
BIT_SIZE 60, 73
BLANK 47-48
blank 6-7, 9, 20, 28-29, 42, 46-47, 49,
52, 59, 67, 69
BLOCK 82-85
BLOCK DATA 11, 43, 92-94
BTEST 60, 90
C
CALL 14, 17, 19-20, 46, 51-53, 64, 70, 83, 89
CASE 9, 10, 19, 50, 94
CEILING 48
CHAR 58
CHARACTER 8, 17, 43, 52, 55
character string 2, 4, 6-8, 27, 47, 49, 52, 56,
58-59, 69, 72, 75, 80
CLOSE 46
CMPLX 5, 48, 58
comment 3, 9, 27, 78, 90
comment character 9, 27, 49, 75, 90
COMMON 6, 29, 43, 55, 83, 89
compatible 2, 9, 57
compiler 26, 41, 71, 74-75
COMPLEX 6, 24, 26, 42-43, 49, 59, 73, 75
conditional GOTO-statement 10, 44
conditional statement 9, 45, 75, 78
Connection Machine 74
CONJG 48
CONTAINS 11-12, 20, 80, 93
continuation character 9, 29
CONTINUE 45, 57, 78-79
control variable 44, 45, 51, 57
Convex 74
COS 7, 48
COSH 48
COUNT 62, 75, 88
CRAY 6, 9, 24, 38, 40, 75
CSHIFT 66-68, 75
CYCLE 9-10, 29, 51
CYCLIC 82, 84-85
D
DATA 4, 27, 44, 55
data flow analysis 27, 72
DATE_AND_TIME 69
DBLE 48
DBOS 73
DEALLOCATE 21, 53, 91
DEC 6, 9, 27, 38-39, 43, 48, 71, 73-74, 88
declaration (specification) 4-6, 11, 13, 16, 21-24, 26,
28-29, 32-35, 37,42-43, 49-50, 53, 55, 72, 75, 80
DEFAULT 10
default argument 14-15, 51, 90
DELETE 48
DELIM 48
DIGITS 60, 73
-----------------------------------------------------------------------
- 102 -
DIM (argument) 36, 37, 87-88, 90
DIM (numeric function) 48
DIMENSION 4, 6, 37, 43, 49, 78-80, 91-92
DIRECT 48
directive 81-82
directory 15, 72
DISTRIBUTE 82-83, 85, 89
DO-loop 2, 9, 51, 57, 75, 78-80, 85-87
DOT_PRODUCT 62, 75
double precision 5, 24, 77
DOUBLE PRECISION 5, 24, 26, 42, 43, 75
DPROD 48
DYNAMIC 84, 90
dynamic memory allocation 10, 23, 25, 37
E
ELSE 15, 42, 55, 79-80
ELSEWHERE 52, 90
END 8, 9, 45, 48-49, 72, 75, 79-80, 86
ENDFILE 46
ENTRY 43
EOR 48
EOSHIFT 67-68, 75
EPSILON 60, 73
EQUIVALENCE 43, 89
ERR 48
error message 57, 72
exercise 3-4, 6-8, 10, 14-15, 17, 20, 23
EXIST 48
EXIT 9-10, 29, 51
EXP 7, 48
EXPONENT 61
extended DO loop 2, 57
extent 19, 52, 92
external function 11, 16, 31-33, 43, 74, 80, 88, 90, 92
external subroutine 11, 52, 92
EXTERNAL 31, 43, 92
EXTRINSIC 90
F
FILE 48
fix form 3, 6, 9, 26, 28-29, 49, 71-72, 75, 82
FLOOR 48
FMT 48
FORALL 81, 85-87, 90
FORM 48
FORMAT 8, 25, 46, 57, 72, 78-79
format directed input/output 46-48
FORMATTED 48
forward compatibility 57
FRACTION 61
free form 3, 6, 9, 27, 29, 49, 71-72, 82
FTN90 71, 73
function 1, 3-5, 11, 15-17, 21, 24, 26, 29,
30, 32, 34-37, 42-44,
46, 49, 52-54, 58-70, 74, 76-80, 90, 92, 95-97
FUNCTION 12-13, 16, 31-34, 43, 49-50, 52, 80
-----------------------------------------------------------------------
- 103 -
G
generic 17, 90, 92
GO (part of GO TO) 44, 49
GO TO 44, 49, 94
GOTO-statement 10, 44-45, 51, 57
GRADE_DOWN 88
GRADE_UP 88
H
Hewlett Packard 48, 71, 74
hexadecimal number 26, 47, 55, 75, 90
high performance Fortran 2, 81-90
Hollerith-constant 2, 47, 57
HPF 2, 39, 57, 81-90
HUGE 26, 60, 73
I
IACHAR 58
IALL 88
IAND 60, 88, 90
IANY 88
IBCLR 61, 90
IBITS 61, 90
IBM 1, 6, 26, 28, 71, 73, 74
IBSET 61, 90
ICHAR 58
IEOR 61, 88, 90
IF-statement 9, 15, 30-31, 42, 45, 57, 78-80, 94
ILEN 88
I-list 15, 72
IMPLICIT 43
IMPLICIT NONE 4, 8, 13-20, 28, 43, 49, 75, 90, 93
IN 11, 30-34, 51, 95-96
intrinsic function 54, 58-69, 81, 87-88, 90, 93
intrinsic subroutine 69-70, 93
INCLUDE 9, 49, 75, 90
INDEPENDENT 86-87
INDEX 26, 59
INOUT 11, 18, 95-96
input 25, 46, 47, 88
input variable 11, 70
INQUIRE 46, 75
INT 5, 48, 58
INTEGER 5, 24, 25, 43, 50
INTENT 11, 15, 17-18, 32-33, 51, 92-93, 95-96
interface 1, 14-15, 35, 54, 70, 90-91, 92
INTERFACE 12, 14, 15, 17-18, 31-33, 37, 52, 54, 63, 90-91, 93
internal function 11, 16, 80, 93
internal subroutine 11, 93
INTRINSIC 44, 54, 92-94
IOLENGTH 48
IOR 61, 88, 90
IOSTAT 48
IPARITY 88
ISHFT 61, 90
ISHFTC 61, 90
-----------------------------------------------------------------------
- 104 -
K
KEEP 48
keyword 14, 51, 90
kind 4-6, 24, 27, 58-61, 70, 73
KIND 4, 24, 26-28, 59, 73, 90
L
LBOUND 37, 54, 63
leading dimension 25
LEADZ 88
LEN 59
LEN_TRIM 59
list directed input/output 25
LGE 59
LGT 59
library 26, 27, 90
linker 26, 71-72
LINK77 71-72
LLE 59
LLT 59
LOCAL 90
LOG 48
LOG10 48
LOGICAL 4, 26, 28, 43, 60
logical IF-statement 26
loop 94
M
main program 3, 11, 13, 18-19, 33, 43, 45, 80, 92, 95-97
matrix 6, 20, 38, 52
MASK 26
MATMUL 19, 20, 62, 75
MAX 7, 30, 48
MAXEXPONENT 60, 73
MAXLOC 35, 68, 75, 87, 90
MAXVAL 62, 75, 87-88
MERGE 63-64, 75
MIL-STD-1753 55, 81, 90
MIN 30, 48
MINEXPONENT 60, 73
MINLOC 68, 75, 87, 90
MINVAL 62, 75, 87-88
MOD 48
module 11, 13-15, 18, 40, 55-56, 71-72, 74, 90, 94
MODULE 11-13, 18
MODULO 48
MS-DOS 28, 40, 71, 73
MVBITS 70, 90
-----------------------------------------------------------------------
- 105 -
N
NAG 3, 6, 8, 26-28, 40, 48, 71-73
NAME 48
NAMED 48
NAMELIST 25, 52
NCOPIES 36, 37, 59, 66
NEAREST 61
nested 94
NEW 48
new line 25
NeXT 71
NEXTREC 48
NINT 48, 58
NML 48
non-advancing input/output 25
NONE (in IMPLICIT NONE) 4, 8, 13-20, 25, 28, 31-34, 43, 49, 63, 90
NOT (not .NOT.) 60, 90
NULL 47-48
NULLIFY 21
NUMBER 48
NUMBER_OF_PROCESSORS 88
O
obsolescent 2, 57, 73
octal number 26, 47, 55, 75, 90
OLD 48
ONLY 13
OPEN 46, 47, 75
OPENED 48
OPERATOR 12, 93
optimization 16, 40
OPTIONAL 15, 30-31, 51, 92
OUT 11, 30-34, 96
output 6, 15, 19, 25, 29, 37, 40, 42,
46-48, 55, 72, 76, 78
output variable 11, 16, 20, 70
-----------------------------------------------------------------------
- 106 -
P
Pacific-Sierra 74
PACK 65, 75
PAD 48, 66
parallel 1, 36, 38, 53, 57, 75, 81, 86, 88
parameter (argument) 95-97
PARAMETER 4, 8, 13, 27-28, 43, 49, 92
Parasoft 74
PARITY 88
Pascal 40, 42, 76, 95
PAUSE 45, 57
PC 28, 40, 71, 73
pointer 11, 20-21, 37, 55, 75
POINTER 21, 37, 92
POPCNT 88
POPPAR 88
POSITION 48
precision 24
PRECISION 24, 50, 60, 73
PREFIX 88
PRESENT 15, 30, 58
PRINT 8, 46, 78
PRIVATE 13, 55, 92
PROCEDURE 12
PROCESSORS 82, 84
PROCESSORS_SHAPE 88
PROCESSOR_VIEW 90
PRODUCT 62, 75, 88
PROGRAM 8, 33, 43
program library 25-26
program unit 3, 11, 51, 71, 94
PUBLIC 13, 92
punched card 3, 42, 49, 78
PURE 90
-----------------------------------------------------------------------
- 107 -
Q
QUOTE 48
R
RADIX 60, 73
RANDOM_NUMBER 70
RANDOM_SEED 70
rank 19, 52, 92
RANGE 60, 73
READ 46, 48, 78-80
READWRITE 48
REAL 4, 5, 28-29, 43, 48, 49, 58
REALIGN 85, 89-90
REC 48
RECL 48
recursion 16, 42, 49, 52, 94
RECURSIVE 16, 31-32, 37, 52
REDISTRIBUTE 85, 89-90
REPEAT 59
repetition 27
REPLACE 48
RESHAPE 53, 63-64, 66, 68, 75
RESULT 16, 20, 32-34, 52
result variable 16, 20, 94
RETURN 46
REWIND 46, 48
RRSPACING 61
-----------------------------------------------------------------------
- 108 -
S
SAVE 44, 53, 92
SCALE 26, 61
SCAN 59
SCATTER 88
SCRATCH 48
segmentation error 15, 54
SELECT 29, 50
SELECTED_INT_KIND 5, 59-60
SELECTED_REAL_KIND 4, 24, 28, 50, 72
SET_EXPONENT 62
separator 49
sequence association 89
SEQUENCE 55, 89
SEQUENTIAL 48
shape 19, 36, 52, 63, 66, 67, 75, 92
SHAPE 63
SIGN 48
Silicon Graphics 48, 71, 74
SIN 7, 48
single precision 24, 27
SINH 48
SIZE 26, 34, 48, 63
source code 3, 9, 27, 49
SPACING 62
SPARC 26, 71, 73-74
specification (declaration) 4-6, 11, 13, 16, 21-24, 26, 28-29,
32-35, 37, 42-43, 49-50, 53, 55, 72, 75, 80
SPREAD 36, 37, 54, 66, 75
SQRT 48, 52, 78-80
standard 1, 7, 28, 38, 47, 49, 81
statement function 78-80, 94
statement number variable 44, 57
STATUS 48
storage association 89
STOP 45, 78-79
subprogram 11, 43-46, 91
subroutine 1, 3, 11, 16-17, 42-43, 45-46, 58, 69, 70, 94-97
SUBROUTINE 14-15, 17-20, 30-31, 34-35, 37,
43, 51, 53-54, 63-64, 95-96
SUFFIX 88
SUM 19, 26, 54, 62, 75, 88
Sun 6, 9, 13, 26, 27, 48, 71, 74
Swedish characters 6, 42, 72
syntax error 27, 28
SYSTEM_CLOCK 69
-----------------------------------------------------------------------
- 109 -
T
TAN 7, 48
TANH 48
target 21-22, 69
TARGET 21, 37, 92
TEMPLATE 83-84
TINY 26, 60, 73
THEN 42, 79-80
TO (part of ASSIGN) 44
TO (part of GO TO) 44, 49
TRANSFER 54, 60
transition 25, 74
TRANSPOSE 54, 75, 68
TRIM 26, 52, 59
TYPE 12, 22, 26, 55, 90
U
UBOUND 37, 54, 63
ULTRIX 6, 9, 39, 47, 71, 73-74
unit number 72
UNIX 28, 40, 42, 48, 71-72, 76
UNDEFINED 48
UNFORMATTED 48
UNIT 48, 72
UNKNOWN 48
UNPACK 66, 75
USE 13, 18, 19
V
VAX 9, 71
vector 4, 20, 52, 57, 81
VERIFY 59
W
WHERE 52, 54, 75, 90
WHILE 42, 51, 75, 90
WRITE 8, 29, 37, 46, 48, 79-80
Y
YES 48
Z
ZERO 47-48
-----------------------------------------------------------------------
- 110 -
FORTRAN 90 FOR THE FORTRAN 77 PROGRAMMER
The programming language FORTRAN is the principal language for
scientific and technical computations. It was developed originally in
1954 and has been revised several times, to FORTRAN IV and Fortran 77.
It has recently been carefully revised, resulting in a modern and
powerful language, Fortran 90.
The intent with this new standard is to make Fortran into a useful
and efficient language for the scientific and technical computations
also towards the end of this decade. The new version contains powerful
new features for the treatment of vectors and matrices, several new
possibilities to specify the precision, access to environment
parameters, intrinsic functions for manipulating floating point numbers,
internal procedures and new specifications for storage and interfacing.
In addition there is a new improved layout of the source code, new
conditional statements, recursion and dynamic memory allocation.
Nothing was removed, so Fortran 90 contains the whole of Fortran 77, but
offers both easier programming and improved security.
An important requirement at the introduction of a new language is
nowadays that the language should be efficient also on parallel systems.
This book therefore contains an Appendix on the recent proposed addition
to Fortran, High Performance Fortran.
This book assumes that the reader is experienced in Fortran 77.
Those parts of Fortran 90 that are already in Fortran 77 are not
discussed here. The many programming examples and user exercises
illustrate the programming technique and available commands, and makes
it easy to get started with writing programs using the new facilities in
Fortran 90. All examples have been tested on Sun SPARC (UNIX) and DEC
Station (ULTRIX), and on IBM PC with MS-DOS 5 and 6.2.
------------------------------------------------------------------------
- 111 -
BO EINARSSON (born 1939)
Dr. Techn. Bo Einarsson has more than 25 years of experience as a
Fortran programmer and as a teacher in Fortran. He has been a principal
scientific officer and director of research at the Swedish National
Defence Research Institute in Tumba and Umea. His main area of work
there was computational continuum mechanics. He joined Linkoping
University as director of the university computing center and is now an
associated professor with teaching in numerical methods and Fortran
Programming. He is also systems manager at the Swedish National
Supercomputer Centre, where Fortran is the dominating language on the
Cqray.
He has been active in the standardization of Fortran and floating
point arithmetic, and has founded and been the chairman of the
International Federation for Information Processing (IFIP) working group
on numerical software (IFIP WG 2.5). He has been awarded the IFIP
"Silver Core" and is an honorary professor at the Australian Institute
for Coordinated Research.
YURIJ I SHOKIN (born 1943)
Professor Yurij Shokin has been many years with the Computer Center
and the Institute of Theoretical and Applied Mechanics in Novosibirsk
(during the directorship of Academician N N Yanenko) and was director of
the computer centre of the Soviet Academy of Sciences, Siberian
Division, in Krasnoyarsk. His main area of research is computational
fluid dynamics. He is now director of the Institute of Computational
Technologies of the Russian Academy of Sciences, Siberian Division, in
Novosibirsk.
He is a member of IFIP WG 2.5 and a corresponding member of the
Russian Academy of Sciences. He is the author of many scientific books
in Russian, including
Methods of Differential Approximations, 1979 and 1985
Numerical Simulation of Tsunami Waves, 1983
Methods of Interval analysis, 1986
Numerical Experiments of the Tsunami Problem, 1989
Fortran 90
for the Fortran 77 programmer
Copyright C 1993, Bo Einarsson and Yurij Shokin
------------------------------------------------------------------------
Verantwortlich für den Inhalt: NN
Datum der letzten Änderung: 12.09.1996