2

今、私はモジュールで、シードを使用して関数のループで使用する乱数を生成していますが、その関数を呼び出すたびに同じ乱数が生成されるという問題に直面しています(シードは明らかに同じであるため) )しかし、シリーズを継続する必要があるか、少なくとも呼び出し間で異なる必要があると想定されています。1つの解決策は、メインプログラムがモジュールで使用される新しいシードを提供することですが、別のエレガントな解決策があると思います. 多くの人からの提案により、Mersenne Twister generatorを使用しています。

追加した

モジュール内の私の関数 (関数のパッケージです) は、シードによって生成された乱数を使用してそのような 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 も使用して) メイン プログラムの contains に配置し、contains セクションと 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 でこれを見つけるのに 10 分もかかりませんでした。

于 2013-09-12T19:42:39.620 に答える
0

獲得したポイントを回復するために、私は自分で答えを見つけなければなりませんでした。これが(1時間の試行の後)です

メインプログラムは

    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 コードをモジュールに便利に含めるように変更しました。メインプログラムとは別にメトロポリスのテスト。メイン プログラムは、トライアルが受け入れられたかどうかを知りたいだけです。私の解決策は、メインプログラムをより明確にします]

于 2013-09-17T18:06:20.160 に答える