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 its translation into English is in progress. 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 This version is dated 28 June 1994. -------------------------------------------------------------------- Preface v Introduction 1 1. To observe at 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 25 15. Use of program libraries 26 16. Peculiarities in the language Fortran 90 27 Answers and comments to the user exercises 28 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 71 7. Other Fortran 90 systems 74 8. Fortran 90 properties in Cray CF 77 75 9. The historical development of Fortran 76 10. High Performance Fortran 81 11. Explanations of certain terms 91 12. Argument association 95 13. Index 98 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 constructions, the other is that Fortran 90 gives us a more strict check of the program. This means a greater security in the programs obtained. During the end of 1993 Fortran 90 will be offered by most computer manufacturers and the language will be a success in 1994. 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. To observe at 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 consequent, 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. 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) ---------------------------------------------------------------------- - 6 - 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 constructions 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 construction 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 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 ---------------------------------------------------------------------- - 15 - 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. ---------------------------------------------------------------------- - 25 - 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. 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) ---------------------------------------------------------------------- - 26 - 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. 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. ---------------------------------------------------------------------- - 27 - 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 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. 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. ---------------------------------------------------------------------- - 28 - 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. 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 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. ---------------------------------------------------------------------------- - 30 - (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 ! Initializing M N = SIZE (B) M (1:N, 1:N) = A M (1:N, N+1) = B ---------------------------------------------------------------------- - 35 - ! 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 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. ---------------------------------------------------------------------- - 38 - 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, McGraw-Hill, New-York 1990. ISBN 0-07-000248-7. $34.95. 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. 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. 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 for DEC-programmers. Very expensive. ---------------------------------------------------------------------- - 39 - 9. 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. 10. 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. 11. 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. 12. T. M. R. Elles: 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. 13. 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. 14. 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. 15. 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. ---------------------------------------------------------------------- - 40 - 16. 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. 17. 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. 18. 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. 19. 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. 20. 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. 21. 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. 22. 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 development of Fortran 90. The 1993 printing contains a very complete application example. 23. 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. 24. NAG: NAGWare FTN90 Compiler (MS-DOS), NP2401/2488, 2nd Edition, Revision A, December 1992. ISBN 1-85206-080-8. 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 21). 25. 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. Just published, so I have not read it yet. 26. 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. 27. 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. ---------------------------------------------------------------------- - 41 - 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. ---------------------------------------------------------------------- - 48 - 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). 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 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. ---------------------------------------------------------------------- - 49 - 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