2

现在我面临的问题是，在一个模块中，我使用种子生成随机数以在函数的循环中使用，但每次调用该函数时，都会生成相同的随机数（因为种子显然是相同的) 但假设它必须继续该系列，或者至少它必须在调用之间有所不同。一种解决方案可能是主程序提供了一个新的种子以在模块中使用，但我认为可能还有另一种优雅的解决方案。根据很多人的建议，我正在使用Mersenne Twister发电机。

添加

我在我的模块中的函数（它是一个函数包）使用种子生成的随机数进行这样的 Metropolis 测试，由于某种原因，如果我把

    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

但是，如果我将该函数放在主程序的包含中（也使用 mtmod）并在包含部分和对 myfunc 的调用之前调用 sgrnd(4357)，那么现在一切都可以很好地编译和运行。为了清楚起见，我不想把那么长的函数放在主程序中，它有 70 行代码，但我似乎没有逃脱的余地。请注意，种子被调用一次。模拟现在具有物理意义，但付出了代价。

4

3 回答 3

4

我总是使用这个子程序（我正在运行蒙特卡罗模拟），在主程序的开头调用它，它应该可以完成这项工作：

（来源：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
于 2013-09-13T07:57:34.167 回答
2

您可以在此处找到一个使用系统时间重新播种随机数生成器的子程序。不必每次调用 random_number() 时都执行此操作，只需在每次重新启动程序时执行此操作即可。

老实说，我用 Google 不到十分钟就找到了这个。

于 2013-09-12T19:42:39.620 回答
0

为了恢复我的分数，我不得不找到自己的答案，在这里（经过一小时的尝试）

主程序是

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

这是我的模块

    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

测试我的解决方案并投票给我 ;) [当然，如您所见，我修改了 mt.f90 代码以方便地包含在我的模块中，因此我可以将主程序与随机数生成部分分开，这样我就可以做到一个 Metropolis 测试除了主程序。主程序只想知道试验是否被接受。我的解决方案确实使主要程序更加清晰]

于 2013-09-17T18:06:20.160 回答