(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 ENDREMARK: 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+06Numbers 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 DOThe 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 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 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_QUADThe 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_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_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.
|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 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 :: 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 EXTERNAL statement is no longer required, but you get a nested INTERFACE.
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_EQNSThe 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 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 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.
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 DOmay 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 DOThis 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