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, NAG's Fortran 90, for MS-DOS and UNIX computers. 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) 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 may give 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
          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)

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, but might not be detected until execution time. You can instead write

            WRITE(*,*) ' HI '
or
            WRITE(*,'(A)') ' HI '
which both write out the text HI on the standard unit for output. If you wish to give the text, which you wish to print, directly where the output format is to be given, this can be done with either apostrophe editing as
            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+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 DO
The English word sum is not suited as the variable name in this case, since this is also an intrinsic function. Summa is the Swedish word for sum.

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

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

(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 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_QUAD
The file test_qua.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.
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.
       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. 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_qua.f90 contains the main program and the function f(x), while the second quad.f90 contains the recursive function.

% f90 test_qua.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.
ComputerTimePrecision
secondsdecimal digits
PC Intel 486 SX 25 74.86
PC Intel 486 DX 50 2.756
PC Intel Pentium 20000.266
Sun SPARC SLC 2.506
Sun SPARC station 100.586
Cray Y-MP 0.1914
Cray C 90 0.1314
DEC Alpha 3000/900 0.066
DEC Alpha 3000/900 0.0615
DEC Alpha 3000/900 3.3233
The time on the Cray is not especially small, since this program could not utilize vectorization, which is one of the main advantages with Cray.

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.

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_qu2.f90 contains the main program and the function f(x), while the second quad2.f90 contains the recursive function.

(11.1)

SUBROUTINE SOLVE_SYSTEM_OF_LINEAR_EQNS(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_LOC(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_EQNS
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).

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.

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.

The DO-loop

        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
may cause a signal on out-of-bound index from an empty sum when K = N and index checking (not standard) is on. With the following slight change this potential problem is avoided.
        X (N) = M (N, N+1) / M (N, N)
        DO K = N-1, 1, -1
               X (K) = (M (K, N+1) - &
                 SUM (M (K, K+1:N)* X (K+1:N)) ) / M (K, K)
        END DO
This version is now included in the single precision and double precision routines stored as Fortran files.

(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 matrix in the subroutine. In this way we get a dynamically allocated matrix.
      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: 11 December 2001
boein@nsc.liu.se