calculation of the integral method trapeze

Question

I have already current code, but it still not working. If code is correct, please help how I can compile it. I had tried it to compile so:

gfortran trap.f -fopenmp

PROGRAM TRAP
USE OMP_LIB
DOUBLE PRECISION INTEG, TMPINT  
DOUBLE PRECISION A, B      
PARAMETER (A=3.0, B=7.0)  
INTEGER N      
PARAMETER (N=10)
DOUBLE PRECISION H         
DOUBLE PRECISION X
INTEGER I
DOUBLE PRECISION F         
H = (B-A)/N
INTEG = 0.0
TMPINT = 0.0
!$omp parallel firstprivate(X, TMPINT) shared(INTEG)
!$omp do
DO 10 I=1,N-1,1
X=A+I*H
TMPINT = TMPINT + F(X)
10   CONTINUE
!$omp end do
!$omp critical
INTEG = INTEG + TMPINT
!$omp end critical
!$omp end parallel 
NTEG = (INTEG+(F(A)+F(B))/2.0)*H
PRINT *, "WITH N=", N, "INTEGRAL=", INTEG
END
FUNCTION F(X)
DOUBLE PRECISION X
F = X / (X + 1) * EXP(-X + 2)
END

Compiler gives following problems:

[https://i.stack.imgur.com/QPknv.png][1]

[https://i.stack.imgur.com/GYkmN.png][2]

Answer 1

Your program has a suffix .f, so gfortran assumes that the code is in fixed format and complains that many statements are "unclassifiable". To fix this, change the file name to trap.f90 and compile as gfortran -fopenmp trap.f90 to assume free format. There are also other problems: one is that the return type of function F(X) does not match with the type declared in the main program, so F(X) needs to be modified as

FUNCTION F(X)
implicit none                     !<--- this is always recommended
DOUBLE PRECISION X, F             !<--- add F here
F = X / (X + 1) * EXP(-X + 2)
END

Another issue is that NTEG is probably a typo of INTEG, so it should be modified as

INTEG = (INTEG+(F(A)+F(B))/2.0)*H

(this is automatically detected if we have implicit none in the main program). Now running the code with, e.g. 8 threads, gives

$ OMP_NUM_THREADS=8 ./a.out
WITH N=          10 INTEGRAL=  0.28927708626319770

while the exact result is 0.28598... Increasing the value of N, we can confirm that the agreement becomes better:

WITH N=         100 INTEGRAL=  0.28602065571967972
WITH N=        1000 INTEGRAL=  0.28598803555916535
WITH N=       10000 INTEGRAL=  0.28598770935198736
WITH N=      100000 INTEGRAL=  0.28598770608991503

---

BTW, it is probably easier to use the reduction clause to do the same thing, for example:

INTEG = 0.0

!$omp parallel do reduction(+ : integ) private(x)                                   
DO I = 1, N-1                                                                       
    X = A + I * H
    INTEG = INTEG + F( X )
ENDDO
!$omp end parallel do                                                               

INTEG = (INTEG+(F(A)+F(B))/2.0)*H

Answer 2

Your code is in fixed form (.f). Therefore, you must code by the rules of the fixed format: The first six characters on each line have a special meaning and should be blank unless you specify a comment in the first position, a line continuation (sixth position), or statement labels 10.

If you format your code accordingly, the compiler complains about a mismatch in the return value of F(X). As you do not use implicit none, the type is defined by the first letter of the function, and F maps to a (single precision) real. So you need to specify the return type explicitly.

Then the code looks like:

      PROGRAM TRAP
      USE OMP_LIB
      DOUBLE PRECISION INTEG, TMPINT  
      DOUBLE PRECISION A, B      
      PARAMETER (A=3.0, B=7.0)  
      INTEGER N      
      PARAMETER (N=10)
      DOUBLE PRECISION H         
      DOUBLE PRECISION X
      INTEGER I
      DOUBLE PRECISION F         
      H = (B-A)/N
      INTEG = 0.0
      TMPINT = 0.0
c$omp parallel firstprivate(X, TMPINT) shared(INTEG)
c$omp do
      DO 10 I=1,N-1,1
      X=A+I*H
      TMPINT = TMPINT + F(X)
10    CONTINUE
c$omp end do
c$omp critical
      INTEG = INTEG + TMPINT
c$omp end critical
c$omp end parallel 
      INTEG = (INTEG+(F(A)+F(B))/2.0)*H
      PRINT *, "WITH N=", N, "INTEGRAL=", INTEG
      END

      DOUBLE PRECISION FUNCTION F(X)
      DOUBLE PRECISION X
      F = X / (X + 1) * EXP(-X + 2)
      END

[Please note that I also fixed the NTAG = line into INTEG= as I believe this is intended. I did not check the code for validity. ]