(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
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, 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.
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%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 SOLVE
The 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 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.
| Computer | Time | Precision |
|---|---|---|
| seconds | decimal digits | |
| PC Intel 486 SX 25 | 74.8 | 6 |
| PC Intel 486 DX 50 | 2.75 | 6 |
| PC Intel Pentium 2000 | 0.26 | 6 |
| Sun SPARC SLC | 2.50 | 6 |
| Sun SPARC station 10 | 0.58 | 6 |
| Cray Y-MP | 0.19 | 14 |
| Cray C 90 | 0.13 | 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 :: 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.
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.
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.
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