generating random numbers in a Fortran Module

2019-05-19 08:04发布

问题:

Now I am facing the problem that in a module, with a seed I am generating random numbers to be used in a loop of a function but each time I call that function, the same random numbers are generated (because the seed is obviously the same) but it's supposed that it must continue the series or at least it must be different between calls. One solution could be that the main program gives a new seed to be used in the module but I think it there could be another elegant solution. I am using Mersenne Twister generator by suggestion of many people.

Added

My function in my module (it is a package of functions) escentially makes such a Metropolis test using random numbers generated by a seed, for some reason 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 that function in the contains of the main program (using mtmod too) and call sgrnd(4357) before contains section and the calls to myfunc, now everything compile and run nicely. For clarity, I didn't want to put that long function in the main program, it has 70 lines of code, but it seems I have no escape. Notice that so, the seed is once called. The simulations have now physical meanings but with that price payed.

回答1:

I always used this subroutine (I'm running a MonteCarlo simulation), call it in the beginning of your main program and tis 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


回答2:

You can find here a subroutine that uses system time to re-seed the random number generator. You shouldn't have to do this every time you call random_number(), just each time you re-start the program.

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



回答3:

In order to recover my points taken, I was obliged to find my own answer, here it is (after one hour of tries)

main program is

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

here's 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(a-z)
        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 w-r 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,N-M-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=N-M,N-2
          y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
          mt(kk)=ieor(ieor(mt(kk+(M-N)),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

test my solution and vote me up ;) [of course, as you see, I modified mt.f90 code to be included conveniently in my module, so I can keep separately the main program from the randon numbers generation part, so I can do a Metropolis test aside the main program. The main program just wants to know if a trial was accepted or not. My solution does give more clarity to the main progam]