Quasi Monte Carlo Integration

According to the theory, the error of quasi Monte
Carlo integration is about O(N^-1), contrasting
to O(N^-1/2) for ordinary Monte Carlo method.

It is noted that O(N^-1) has huge speed advantage
comparing to Monte Carlo method.

The following is the fortran 90 code for quasi Monte
Carlo method with Halton sequence in 4D. Notice
here we do not use parallel version for Halton
sequence.

MODULE GLOBE

REAL*8,PARAMETER::LAMBDA=3.0d0
REAL*8,PARAMETER::GAMMA=1.0d0
REAL*8,PARAMETER::PI=3.141592653589793238462643383279592884197169399375d0
REAL*8,PARAMETER::PI2=PI*2.0d0

ENDMODULE GLOBE



Program quasimontecarlo
use GLOBE
external F
include 'mpif.h'
REAL*8:: a,b,c(1:4),func,kx,ky,qx,qy
INTEGER*8 ,DIMENSION(1:4) :: N1,NRA1, NRB1
REAL*8,DIMENSION(1:4) :: RS1
REAL*8::RESULTI,x,y,z,G,REALvalue
integer:: myid,numprocs,ierr1,timer0,timer
INTEGER*8::ISTEP,i,j
integer::L
call MPI_INIT( ierr1 )
call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr1 )
call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr1)
print*,'Lambda=',Lambda, ' Gamma=',Gamma
L=myid+1
timer0=mpi_wtime()
NRA1(:)=0
NRB1(:)=1
N1(1)=2
N1(2)=3
N1(3)=5
N1(4)=7
RESULTI=0.0d0
ISTEP=50000000000
! print*,'test 2'
DO j=1,20
CALL HALTON( N1,NRA1, NRB1,RS1)
END DO
! print *, 'test 1'

DO j=1,ISTEP
CALL HALTON( N1,NRA1, NRB1,RS1)
kx=RS1(1)*PI2
ky=RS1(2)*PI2
qx=RS1(3)*PI2
qy=RS1(4)*PI2
! print*,x,y,z
RESULTI=RESULTI+F(kx,ky,qx,qy,L)
END DO
RESULTI= RESULTI/dble(ISTEP)
timer=mpi_wtime()
print*, myid+1, RESULTI
timer=mpi_wtime()
print*, myid+1, RESULTI

! print*, 'FOR single ISTEP=',ISTEP, (timer-timer0), 'second with result ',RESULTI

call mpi_finalize(ierr1)

END PROGRAM quasimontecarlo




FUNCTION F(kx,ky,qx,qy,L) RESULT(F_RESULT)
USE GLOBE
Integer:: L
REAL*8::x,y,z,F_RESULT,alpha,f1,f2,kx,ky,qx,qy,g1,g2,o1,o2,mk
REAL*8:: fejer1,fejer2,dl ,dc

! p1=(dsin(x)**2)*(dsin(y)**2)*(dsin(z)**2)
! p2= (dcos(x)**2)*(dcos(y)**2)*(dcos(z)**2)
! The function ATAN2(Y,X) = arctan(y,x) will return a value in (-pi, pi].
! If Y is positive the result will be positive. If Y is zero the result will
! be zero if X is positive, and pi if X is negative. If Y is negative the
! result will be negative. If X is zero the result will be plus or minus pi/2.
! Both X and Y are not permitted to be zero simultaneously. The purpose of
! the function is to avoid division by zero.
dl=dble(L)
f1=LAMBDA-(dcos(kx)+dcos(ky))
f2=LAMBDA-(dcos(kx+qx)+dcos(ky+qy))
g1=GAMMA* (dsin(kx)+dsin(ky))
g2=GAMMA* (dsin(kx+qx)+dsin(ky+qy))
o1=Dsqrt(f1*f1+g1*g1)
o2=Dsqrt(f2*f2+g2*g2)
If (( o1 .eq. 0.0d0 ) .OR. (o2 .eq. 0.0d0)) THEN
mk=1.0d0
else
mk=1.0d0- (f1*f2+g1*g2)/(o1*o2)
end if
dc=Dcos(qx)
If ( dc .eq. 1.0d0) THEN
fejer1=dl*dl
else
fejer1=(Dcos(qx*dl)-1.0d0)/(dc-1.0d0)
end if
dc=Dcos(qy)
If ( dc .eq. 1.0d0) THEN
fejer2=dl*dl
else
fejer2=(Dcos(qy*dl)-1.0d0)/(dc-1.0d0)
end if


F_RESULT=fejer1*fejer2*mk
F_RESULT=fejer1*fejer2*mk

ENDFUNCTION F




!----------------Halton sequence-----------------------
! The initialization is:
!NRA1 = 0
!NRB1 = 1
!N1 = p, the base one wishes to use. the pth prime
!DIM is the dimension
SUBROUTINE HALTON( N1,NRA1, NRB1,RS1)
INTEGER :: i
INTEGER*8 :: NXA1,NXB1
INTEGER*8,DIMENSION(1:4) :: N1,NRA1, NRB1
REAL*8,DIMENSION(1:4) :: RS1
DO i=1,4
NXA1 = NRB1(i) - NRA1(i)
! Special case:
IF (NXA1 .EQ. 1) THEN
NRA1(i) = 1
NRB1(i) = NRB1(i)*N1(i)
GOTO 19
END IF
! General case:
NXB1 = NRB1(i)/N1(i)
10 IF (NXA1 .LE. NXB1) THEN
NXB1 = NXB1/N1(i)
GOTO 10
ELSE
NRA1(i) = (1 + N1(i))*NXB1- NXA1
END IF
19 RS1(i) = dble(NRA1(i))/dble(NRB1(i))
END DO
RETURN
END SUBROUTINE HALTON
!--------------------------------------------