在Fortran模块中生成随机数 [英] generating random numbers in a Fortran Module
问题描述
现在我面临的问题是，在一个模块中，使用种子生成随机数用于函数循环，但每次调用该函数时，都会使用相同的随机数（因为种子显然是相同的），但是它假设它必须继续这个系列，或者至少它在调用之间必须是不同的。一种解决方案可能是主程序为模块提供了一个新的种子，但我认为这可能是另一个优雅的解决方案。
我使用 Mersenne Twister 生成器，由许多人建议。
新增
我的模块中的函数（它是一个函数包）基本上使用由种子生成的随机数进行Metropolis测试，出于某种原因，如果我放入
module mymod
使用mtmod
调用sgrnd（4357）！< - 此行导致编译错误
包含
myfunc（args）
隐式无
//声明等
！call sgrnd（4357）< - 如果我把这个调用放在这里compilator说ok，
！但是每次调用这个函数时重新开始随机数序列:(
....
！下面的部分在循环中
if（prob< grnd（） ）然后
！grnd（）是随机数生成的
返回
else继续测试到循环结束
结束myfunc
但是如果我把那个fu放了在包含主程序（使用mtmod）中调用sgrnd（4357），然后调用myfunc，现在所有的东西都编译并运行得很好。为了清楚起见，我不想在主程序中使用这个长函数，它有70行代码，但似乎我没有逃脱。请注意，种子曾经被调用过。这些模拟现在具有物理意义，但是支付了这个价格。
为了恢复我的观点，我不得不寻找我自己的答案，在这里（尝试一小时后）
主程序是
$ pre> （$，*）x + writerandnum（）
write（*，*）x + writerandnum（）
write（*，*）x + writerandnum（）
结束程序callrtmod
这里是我的模块
模块mymod
隐式无
！------------- mt变量-------------
！默认种子
整数，参数:: defaultsd = 4357
！周期参数
整数，参数:: N = 624，N1 = N + 1
！状态向量的数组
整数，保存，维（0：N-1）:: mt
整数，保存:: mti = N1
！-------- ------------------------------
包含
函数writerandnum
隐式none
real（8）:: writerandnum
writerandnum = grnd（）
！如果您愿意，也可以在这里执行Metropolis测试
end function writerandnum
！初始化子程序
子程序sgrnd（种子）
隐式无
整数，意图（in）:: seed
mt（0）= iand（seed，-1）
do mti = 1，N-1
mt（mti）= iand（69069 * mt（mti-1）， - 1 ）
enddo
！
返回
结束子程序sgrnd
！--------------------------------- ------------------------------------------
！函数grnd在这里
！------------------------------------------- --------------------------------
子程序mtsavef（fname，forma ）
字符（*），intent（in）:: fname
字符，intent（in）:: forma
选择大小写（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
结束选择
关闭（10）
返回
结束子程序mtsavef
子程序mtsaveu（unum，forma）
整数，意图（in）:: unum
字符，意图（in）:: forma
选择案例（forma）
案例（'u'，'U'）
写入（unum）mti
写入（unum）mt
案例默认
写入（unum，*）mti
写入（unum，*）mt
结束选择
返回
结束子程序mtsaveu
子程序mtgetf（fname，forma）
字符（*），intent（in）:: fname
字符，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
结束选择
关闭（10）
返回
结束子程序mtgetf
子程序mtgetu（unum，forma）
整数，意图（in）:: unum
字符，意图（in）:: forma
选择案例（forma）
案例（'u'，'U'）
读取（unum）mti
读取（unum）mt
案例默认
读取（unum，*）mti
读取（unum，*）mt
结束选择
返回
结束子程序mtgetu
！======================================= ========
！随机数字发生器
！ real（8）函数grnd（）
函数grnd！agregue yo
隐式整数（a-z）
real（8）grnd！agregue yo
！周期参数
整数，参数:: M = 397，MATA = -1727483681
！常量向量a
整数，参数:: LMASK = 2147483647
！最低有效位r
整数，参数:: UMASK = -LMASK - 1
！最重要的w-r位
！回火参数
整数，参数:: TMASKB = -1658038656，TMASKC = -272236544
尺寸mag01（0：1）
数据mag01 / 0，MATA /
保存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
！一次生成N个单词
if（mti.eq.N + 1）then
！如果sgrnd（）没有被调用，
调用sgrnd（defaultsd）
！使用默认的初始种子
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 ），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.0 d0）
else
grnd = dble（y）/（2.0d0 ** 32-1.0d0）
endif
返回
结束函数grnd
结束模块mymod
测试我的解决方案并投票给我; [当然，正如你所看到的，我修改了mt.f90代码以方便地包含在我的模块中，因此我可以将主程序与randon数字分开因此我可以在主节目旁边做Metropolis测试。主程序只是想知道是否接受了审判。我的解决方案确实为主程序提供了更多的清晰性。
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.
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]
这篇关于在Fortran模块中生成随机数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
