(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) If 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. Also note that with fix form text in positions 73 to 80 was considered to be a comment.

(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
carriage returns as below.

LOGICAL L L = .FALSE. IF (L) THEN Z = 1.0 ELSE Y=Z END IF ENDREMARK: Also certain Fortran 77 compilers give a warning about the variable

(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`

INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,99)

(2.5) `REAL (KIND=DP) :: E, PI`

REAL (KIND=DP), PARAMETER :: 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.

(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 time but at execution time. You can instead write

WRITE(*,*) ' HI 'or

WRITE(*,'(A)') ' HI 'which both write out the text

WRITE(*, "(' HI ')")or with the obsolescent Hollerith editing

WRITE(*, "(4H HI )")

(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+06Numbers from

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)) ENDThe English word

(7.1) Use the functions `MIN ` and `MAX `
to find the smallest and largest values of all the combinations

A%LOWER * B%LOWER, A%LOWER * B%UPPER, A%UPPER * B%LOWER, A%UPPER * B%UPPERat 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.

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 SOLVEThe very simple integral calculation above can be replaced by the adaptive quadrature in exercise (9.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

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 TRIBONACCIThe 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 ENDand gives the result

(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 Linköping, 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.

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_QUADThe file

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_QUADWe 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.

FUNCTION F(X) RESULT (FUNCTION_VALUE) IMPLICIT NONE REAL, INTENT(IN) :: X REAL :: FUNCTION_VALUE FUNCTION_VALUE = EXP(-X**2) END FUNCTION FNow it is time to do the test on the Sun computer. 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.

If you wish to test the program yourself the source code is directly available
in two files. The first
`test_quad.f90 `
contains the main program and the function *f(x)*, while the second
`quad.f90 `
contains the recursive function.

% f90 test_quad.f90 quad.f90 test_quad.f90: quad.f90: % a.out The integral is 1.7733453512 with error estimate 0.0049186843 with real error 0.0008914471 The integral is 1.7724548578 with error estimate 0.0003375171 with real error 0.0000009537 The integral is 1.7724541426 with error estimate 0.0000356939 with real error 0.0000002384 The integral is 1.7724540234 with error estimate 0.0000046571 with real error 0.0000001192 The integral is 1.7724539042 with error estimate 0.0000004876 with real error 0.0000000000 %I have run this program on a number of different systems and obtained the following timings. If the run is repeated a slightly different timing may be achieved.

Computer | Time | Precision |
---|---|---|

seconds | decimal digits | |

PC Intel 486 SX 25 | 74.8 | 6 |

PC Intel 486 DX 50 | 2.75 | 6 |

Sun SPARC SLC | 2.50 | 6 |

Sun SPARC station 10 | 0.58 | 6 |

Cray Y-MP | 0.19 | 14 |

DEC Alpha 3000/900 | 0.06 | 6 |

DEC Alpha 3000/900 | 0.06 | 15 |

DEC Alpha 3000/900 | 3.32 | 33 |

In the specification above of the `RECURSIVE FUNCTION ADAPTIVE_QUAD `
you may replace the line

REAL, EXTERNAL :: Fwith 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 INTERFACEWith this method an explicit

**Remark.**

The program above was written to illustrate the use of recursive functions and
adaptive techniques, and was therefore not optimized. The main problem is that
the function *f(x)* is evaluated three (or even four) times at each call,
once for each of the present boundary points and twice for the middle point.
Please note that the function values at each of the boundary points were
evaluated already in the previous step.

Thus the obvious change is to include the boundary function values in the list of arguments, and to evaluate the middle point function value only once. In this way the execution time is reduced by a factor of about three.

The revised program is also directly available
in two files. The first
`test_quad2.f90 `
contains the main program and the function *f(x)*, while the second
`quad2.f90 `
contains the recursive function.

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 ! 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_EQUATIONSThe input matrix

Please note that the intrinsic function

Also note that the numbering starts with 1, in spite of that we are looking at the vector with the elements running from

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

At the pivoting process we use the vector

In this first version we only use array sections of vector type at the calculations, but we will now introduce the function

The function

If the variable

(/ 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 4I 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 DOwith

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

Unfortunately it is rather difficult to get the parameters to the intrinsic function

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

(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

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 matrix A. ! Its dimensions were determined in the subroutine, ! the number of elements is available as SIZE(A), ! the extent in each direction as SIZE(A,1) and ! as SIZE(A,2). ! END PROGRAM MAIN_PROGRAM SUBROUTINE SUB(B) REAL, DIMENSION (:,:), POINTER :: B INTEGER M, N ! Here we can assign values to M and N, for example ! through an input statement. ! When M and N have been assigned we can allocate B ! as a matrix. ALLOCATE (B(M,N)) ! Now we can use the matrix B. END SUBROUTINE SUB

Last modified: 6 June 1996