program as05b integer,parameter :: n=50 integer i1,i2,i3,i4,i5 real(8) x1,x2,x3,x4,x5,s,h h=0.07 !Rough estimate of the optimal integration step (based on the formula for the total error in 1D case) s=0.0_8 !Print the limits of integration for x5 and the corresponding function values x5=(1-(n+1)/2.0_8)*h write(*,*) 'Integration limits:' write(*,*) 'x5min=',x5,' f(0,0,0,0,x5min)=',atan(sqrt(1.0_8+5*x5*x5))*exp(-(x5*x5)) x5=(n-(n+1)/2.0_8)*h write(*,*) 'x5max=',x5,' f(0,0,0,0,x5max)=',atan(sqrt(1.0_8+5*x5*x5))*exp(-(x5*x5)) !Compute the integral !$omp parallel do private(x1,x2,x3,x4,x5,i1,i2,i3,i4,i5) reduction(+:s) shared(h) do i1=1,n x1=(i1-(n+1)/2.0_8)*h do i2=1,n x2=(i2-(n+1)/2.0_8)*h do i3=1,n x3=(i3-(n+1)/2.0_8)*h do i4=1,n x4=(i4-(n+1)/2.0_8)*h do i5=1,n x5=(i5-(n+1)/2.0_8)*h s=s+atan(sqrt(1.0_8+x1*x1+2*x2*x2+3*x3*x3+4*x4*x4+5*x5*x5))*exp(-(x1*x1+x2*x2+x3*x3+x4*x4+x5*x5)) enddo enddo enddo enddo enddo !$omp end parallel do write(*,*) 'Result: Integral=',s*h*h*h*h*h end program as05b