Random Number Generation in Fortran Module

Question

Now I am faced with the problem that in the module , with the seed, I generate random numbers that will be used in the function loop, but every time I call this function, the same random numbers (since the seed is obviously the same ), but it is assumed that it should continue the series or, at least, should be different between calls. One solution may be that the main program gives a new seed to be used in the module, but I think it could be another elegant solution. I use the Mersenne Twister generator as suggested by many people.

Added

My function in my module (this is a package of functions) allows me to do such a Metropolis test using random numbers generated by a seed, for some reason the compilation complains if I put

module mymod uses mtmod call sgrnd(4357)!<-- this line causes compilation error contains myfunc(args) implicit none // declarations etc !call sgrnd(4357) <-- if I put this call here compilator says ok, !but re-start random number series each time this function is called :( .... !the following part is inside a loop if (prob < grnd()) then !grnd() is random number generated return else continue testing to the end of the loop cycle end myfunc 

But if I put this function in the main program (using mtmod), and the sgrnd (4357) call before containing the section and myfunc calls, now everything compiles and starts beautifully. For clarity, I did not want to put this long function in the main program, it has 70 lines of code, but it seems that I could not escape. Notice that so, the seed is called. Modeling now has physical meanings, but with this price.

Answer 1

To restore my points, I was obliged to find my own answer, here it is (after one hour of trying)

Main program

  program callrtmod use mymod implicit none real::x x=1.0 write(*,*) x+writerandnum() write(*,*) x+writerandnum() write(*,*) x+writerandnum() end program callrtmod 

here is my module

  module mymod implicit none !-------------mt variables------------- ! Default seed integer, parameter :: defaultsd = 4357 ! Period parameters integer, parameter :: N = 624, N1 = N + 1 ! the array for the state vector integer, save, dimension(0:N-1) :: mt integer, save :: mti = N1 !-------------------------------------- contains function writerandnum implicit none real(8)::writerandnum writerandnum = grnd() !if you please, you could perform a Metropolis test here too end function writerandnum !Initialization subroutine subroutine sgrnd(seed) implicit none integer, intent(in) :: seed mt(0) = iand(seed,-1) do mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) enddo ! return end subroutine sgrnd !--------------------------------------------------------------------------- !the function grnd was here !--------------------------------------------------------------------------- subroutine mtsavef( fname, forma ) character(*), intent(in) :: fname character, intent(in) :: forma select case (forma) case('u','U') open(unit=10,file=trim(fname),status='UNKNOWN',form='UNFORMATTED', & position='APPEND') write(10)mti write(10)mt case default open(unit=10,file=trim(fname),status='UNKNOWN',form='FORMATTED', & position='APPEND') write(10,*)mti write(10,*)mt end select close(10) return end subroutine mtsavef subroutine mtsaveu( unum, forma ) integer, intent(in) :: unum character, intent(in) :: forma select case (forma) case('u','U') write(unum)mti write(unum)mt case default write(unum,*)mti write(unum,*)mt end select return end subroutine mtsaveu subroutine mtgetf( fname, forma ) character(*), intent(in) :: fname character, intent(in) :: forma select case (forma) case('u','U') open(unit=10,file=trim(fname),status='OLD',form='UNFORMATTED') read(10)mti read(10)mt case default open(unit=10,file=trim(fname),status='OLD',form='FORMATTED') read(10,*)mti read(10,*)mt end select close(10) return end subroutine mtgetf subroutine mtgetu( unum, forma ) integer, intent(in) :: unum character, intent(in) :: forma select case (forma) case('u','U') read(unum)mti read(unum)mt case default read(unum,*)mti read(unum,*)mt end select return end subroutine mtgetu !=============================================== !Random number generator ! real(8) function grnd() function grnd !agregue yo implicit integer(az) real(8) grnd !agregue yo ! Period parameters integer, parameter :: M = 397, MATA = -1727483681 ! constant vector a integer, parameter :: LMASK = 2147483647 ! least significant r bits integer, parameter :: UMASK = -LMASK - 1 ! most significant wr bits ! Tempering parameters integer, parameter :: TMASKB= -1658038656, TMASKC= -272236544 dimension mag01(0:1) data mag01/0, MATA/ save mag01 ! mag01(x) = x * MATA for x=0,1 TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) if(mti.ge.N) then ! generate N words at one time if(mti.eq.N+1) then ! if sgrnd() has not been called, call sgrnd( defaultsd ) ! a default initial seed is used endif do kk=0,NM-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) enddo do kk=NM,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(MN)),ishft(y,-1)),mag01(iand(y,1))) enddo y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif y=mt(mti) mti = mti + 1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) if(y .lt. 0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif return end function grnd end module mymod 

check my decision and vote for me;) [of course, as you can see, I modified the mt.f90 code, which should be conveniently included in my module, so I can separately store the main program from the generation part of randon numbers, so I can do a test Metropolis towards the main program. The main program just wants to know if the decision was made or not. My decision really gives more clarity to the main program]

Answer 2

I always used this routine (I run the MonteCarlo simulation), call it at the beginning of your main program, and this should do the job:

(Source: Gfortran 4.6.1 )

 c initialize a random seed from the system clock at every run (fortran 95 code) subroutine init_random_seed() INTEGER :: i, n, clock INTEGER, DIMENSION(:), ALLOCATABLE :: seed CALL RANDOM_SEED(size = n) ALLOCATE(seed(n)) CALL SYSTEM_CLOCK(COUNT=clock) seed = clock + 37 * (/ (i - 1, i = 1, n) /) CALL RANDOM_SEED(PUT = seed) DEALLOCATE(seed) end 

Answer 3

You can find here a routine that uses the system time to reseed the random number generator. You should not do this every time you call random_number (), every time you re-run the program.

Honestly, it didn't take me more than ten minutes to find this with Google.