Naive implementation of parallel Mersenne Twister (MT) pseudorandom number generator
in [Fortran]
|
Prev: Bug in gfortran when reading end of file?
Next: Naive implementation of parallel Mersenne Twister (MT) pseudorandom number generator
From: Gordon Sande on 24 Apr 2008 11:05 On 2008-04-24 11:45:55 -0300, mjm2114(a)columbia.edu said: > Hi there, > I have a question on a naive implementation of a parallel MT that I've > done using the fortran version of MT19937ar.f posted in Prof. > Matsumoto's website. (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ > VERSIONS/FORTRAN/mt19937ar.f) > > First, I setup an MT in the master node and seed it with a single > integer using subroutine init_genrand(seed). I then use the first 624 > outputs from the master node to seed the MT in the first worker node > using subroutine init_by_array(array_of_seeds,624). I continue in the > same manner for all subsequent worker nodes, taking the next 624 > outputs from the master node and using them to seed the MT in all > worker nodes. Once this is done, I have a different MT ready for use > in all the worker nodes. Do you think this is a good approach? I know > that the all blocks of 624 seeds used to set up the worker node MTs > have some correlation (since they are adjacent blocks produced by the > same MT), but given the equidistribution properties and the gigantic > period of the algorithm, wouldn't these correlations be insignificant > for all practical purposes? Would it be better if I used a different > PRNG algorithm for the master node? If that is so, which one would you > recommend? Are gfortran's RAND or RANDOM_NUMBER up to the task? Is there some reason why the second value from the sequence on the first node is not the same as the first value from the second node? (Or some technical variant of a question like this.) Without looking at the details it would seem that the various nodes are just delayed varsions of lower numbered nodes. You would want to have the various nodes be at "independent" places around the full cycle. Having some other PRNG do the node seeding would seem to be more likely to achieve this. > I know I can't expect perfectly uncorrelated streams of outputs in > each node as in serious parallel PRNGs (i.e. SPRNG), but I'm only > going to use the code for small workstations/clusters of 8-16 nodes > (using mpi). > > I also have a question on the output of the MT. The output is supposed > to be signed 32-bit integers, there should be on problem in using > negative integers as seeds for the initialization subroutine right? > > Here's the code > > c > c A C-program for MT19937, with initialization improved 2002/1/26. > c Coded by Takuji Nishimura and Makoto Matsumoto. > c > c Before using, initialize the state by using init_genrand(seed) > c or init_by_array(init_key, key_length). > c > c Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, > c All rights reserved. > c Copyright (C) 2005, Mutsuo Saito, > c All rights reserved. > c > c Redistribution and use in source and binary forms, with or without > c modification, are permitted provided that the following conditions > c are met: > c > c 1. Redistributions of source code must retain the above copyright > c notice, this list of conditions and the following disclaimer. > c > c 2. Redistributions in binary form must reproduce the above > copyright > c notice, this list of conditions and the following disclaimer > in the > c documentation and/or other materials provided with the > distribution. > c > c 3. The names of its contributors may not be used to endorse or > promote > c products derived from this software without specific prior > written > c permission. > c > c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS > c "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT > c LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS > FOR > c A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE > COPYRIGHT OWNER OR > c CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, > SPECIAL, > c EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, > c PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR > c PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY > OF > c LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT > (INCLUDING > c NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS > c SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. > c > c > c Any feedback is very welcome. > c http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html > c email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) > c > c----------------------------------------------------------------------- > c FORTRAN77 translation by Tsuyoshi TADA. (2005/12/19) > c > c ---------- initialize routines ---------- > c subroutine init_genrand(seed): initialize with a seed > c subroutine init_by_array(init_key,key_length): initialize by an > array > c > c ---------- generate functions ---------- > c integer function genrand_int32(): signed 32-bit integer > c integer function genrand_int31(): unsigned 31-bit integer > c double precision function genrand_real1(): [0,1] with 32-bit > resolution > c double precision function genrand_real2(): [0,1) with 32-bit > resolution > c double precision function genrand_real3(): (0,1) with 32-bit > resolution > c double precision function genrand_res53(): (0,1) with 53-bit > resolution > c > c This program uses the following non-standard intrinsics. > c ishft(i,n): If n>0, shifts bits in i by n positions to left. > c If n<0, shifts bits in i by n positions to right. > c iand (i,j): Performs logical AND on corresponding bits of i and > j. > c ior (i,j): Performs inclusive OR on corresponding bits of i and > j. > c ieor (i,j): Performs exclusive OR on corresponding bits of i and > j. > c > c----------------------------------------------------------------------- > c initialize mt(0:N-1) with a seed > c----------------------------------------------------------------------- > subroutine init_genrand(s) > integer s > integer N > integer DONE > integer ALLBIT_MASK > parameter (N=624) > parameter (DONE=123456789) > integer mti,initialized > integer mt(0:N-1) > common /mt_state1/ mti,initialized > common /mt_state2/ mt > common /mt_mask1/ ALLBIT_MASK > c > call mt_initln > mt(0)=iand(s,ALLBIT_MASK) > do 100 mti=1,N-1 > mt(mti)=1812433253* > & ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti > mt(mti)=iand(mt(mti),ALLBIT_MASK) > 100 continue > initialized=DONE > c > return > end > c----------------------------------------------------------------------- > c initialize by an array with array-length > c init_key is the array for initializing keys > c key_length is its length > c----------------------------------------------------------------------- > subroutine init_by_array(init_key,key_length) > integer init_key(0:*) > integer key_length > integer N > integer ALLBIT_MASK > integer TOPBIT_MASK > parameter (N=624) > integer i,j,k > integer mt(0:N-1) > common /mt_state2/ mt > common /mt_mask1/ ALLBIT_MASK > common /mt_mask2/ TOPBIT_MASK > c > call init_genrand(19650218) > i=1 > j=0 > do 100 k=max(N,key_length),1,-1 > mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1664525) > & +init_key(j)+j > mt(i)=iand(mt(i),ALLBIT_MASK) > i=i+1 > j=j+1 > if(i.ge.N)then > mt(0)=mt(N-1) > i=1 > endif > if(j.ge.key_length)then > j=0 > endif > 100 continue > do 200 k=N-1,1,-1 > mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1566083941)- > i > mt(i)=iand(mt(i),ALLBIT_MASK) > i=i+1 > if(i.ge.N)then > mt(0)=mt(N-1) > i=1 > endif > 200 continue > mt(0)=TOPBIT_MASK > c > return > end > c----------------------------------------------------------------------- > c generates a random number on [0,0xffffffff]-interval > c----------------------------------------------------------------------- > function genrand_int32() > integer genrand_int32 > integer N,M > integer DONE > integer UPPER_MASK,LOWER_MASK,MATRIX_A > integer T1_MASK,T2_MASK > parameter (N=624) > parameter (M=397) > parameter (DONE=123456789) > integer mti,initialized > integer mt(0:N-1) > integer y,kk > integer mag01(0:1) > common /mt_state1/ mti,initialized > common /mt_state2/ mt > common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK > common /mt_mag01/ mag01 > c > if(initialized.ne.DONE)then > call init_genrand(21641) > endif > c > if(mti.ge.N)then > do 100 kk=0,N-M-1 > y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK)) > mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) > 100 continue > do 200 kk=N-M,N-1-1 > y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK)) > mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) > 200 continue > y=ior(iand(mt(N-1),UPPER_MASK),iand(mt(0),LOWER_MASK)) > mt(kk)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) > mti=0 > endif > c > y=mt(mti) > mti=mti+1 > c > y=ieor(y,ishft(y,-11)) > y=ieor(y,iand(ishft(y,7),T1_MASK)) > y=ieor(y,iand(ishft(y,15),T2_MASK)) > y=ieor(y,ishft(y,-18)) > c > genrand_int32=y > return > end > c----------------------------------------------------------------------- > c generates a random number on [0,0x7fffffff]-interval > c----------------------------------------------------------------------- > function genrand_int31() > integer genrand_int31 > integer genrand_int32 > genrand_int31=int(ishft(genrand_int32(),-1)) > return > end > c----------------------------------------------------------------------- > c generates a random number on [0,1]-real-interval > c----------------------------------------------------------------------- > function genrand_real1() > double precision genrand_real1,r > integer genrand_int32 > r=dble(genrand_int32()) > if(r.lt.0.d0)r=r+2.d0**32 > genrand_real1=r/4294967295.d0 > return > end > c----------------------------------------------------------------------- > c generates a random number on [0,1)-real-interval > c----------------------------------------------------------------------- > function genrand_real2() > double precision genrand_real2,r > integer genrand_int32 > r=dble(genrand_int32()) > if(r.lt.0.d0)r=r+2.d0**32 > genrand_real2=r/4294967296.d0 > return > end > c----------------------------------------------------------------------- > c generates a random number on (0,1)-real-interval > c----------------------------------------------------------------------- > function genrand_real3() > double precision genrand_real3,r > integer genrand_int32 > r=dble(genrand_int32()) > if(r.lt.0.d0)r=r+2.d0**32 > genrand_real3=(r+0.5d0)/4294967296.d0 > return > end > c----------------------------------------------------------------------- > c generates a random number on [0,1) with 53-bit resolution > c----------------------------------------------------------------------- > function genrand_res53() > double precision genrand_res53 > integer genrand_int32 > double precision a,b > a=dble(ishft(genrand_int32(),-5)) > b=dble(ishft(genrand_int32(),-6)) > if(a.lt.0.d0)a=a+2.d0**32 > if(b.lt.0.d0)b=b+2.d0**32 > genrand_res53=(a*67108864.d0+b)/9007199254740992.d0 > return > end > c----------------------------------------------------------------------- > c initialize large number (over 32-bit constant number) > c----------------------------------------------------------------------- > subroutine mt_initln > integer ALLBIT_MASK > integer TOPBIT_MASK > integer UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK > integer mag01(0:1) > common /mt_mask1/ ALLBIT_MASK > common /mt_mask2/ TOPBIT_MASK > common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK > common /mt_mag01/ mag01 > CC TOPBIT_MASK = Z'80000000' > CC ALLBIT_MASK = Z'ffffffff' > CC UPPER_MASK = Z'80000000' > CC LOWER_MASK = Z'7fffffff' > CC MATRIX_A = Z'9908b0df' > CC T1_MASK = Z'9d2c5680' > CC T2_MASK = Z'efc60000' > TOPBIT_MASK=1073741824 > TOPBIT_MASK=ishft(TOPBIT_MASK,1) > ALLBIT_MASK=2147483647 > ALLBIT_MASK=ior(ALLBIT_MASK,TOPBIT_MASK) > UPPER_MASK=TOPBIT_MASK > LOWER_MASK=2147483647 > MATRIX_A=419999967 > MATRIX_A=ior(MATRIX_A,TOPBIT_MASK) > T1_MASK=489444992 > T1_MASK=ior(T1_MASK,TOPBIT_MASK) > T2_MASK=1875247104 > T2_MASK=ior(T2_MASK,TOPBIT_MASK) > mag01(0)=0 > mag01(1)=MATRIX_A > return > end > > Any advice would be sincerely appreciated. Keep up the good work! > > Manuel
From: Gerry Ford on 24 Apr 2008 18:41
<mjm2114(a)columbia.edu> wrote in message news:6c9d2875-8019-42d0-9389-40850ad575a7(a)x35g2000hsb.googlegroups.com... > Hi there, > I have a question on a naive implementation of a parallel MT that I've > done using the fortran version of MT19937ar.f posted in Prof. > Matsumoto's website. (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ > VERSIONS/FORTRAN/mt19937ar.f) > > First, I setup an MT in the master node and seed it with a single > integer using subroutine init_genrand(seed). I then use the first 624 > outputs from the master node to seed the MT in the first worker node > using subroutine init_by_array(array_of_seeds,624). I continue in the > same manner for all subsequent worker nodes, taking the next 624 > outputs from the master node and using them to seed the MT in all > worker nodes. Once this is done, I have a different MT ready for use > in all the worker nodes. Do you think this is a good approach? I know > that the all blocks of 624 seeds used to set up the worker node MTs > have some correlation (since they are adjacent blocks produced by the > same MT), but given the equidistribution properties and the gigantic > period of the algorithm, wouldn't these correlations be insignificant > for all practical purposes? Would it be better if I used a different > PRNG algorithm for the master node? If that is so, which one would you > recommend? Are gfortran's RAND or RANDOM_NUMBER up to the task? > > I know I can't expect perfectly uncorrelated streams of outputs in > each node as in serious parallel PRNGs (i.e. SPRNG), but I'm only > going to use the code for small workstations/clusters of 8-16 nodes > (using mpi). > > I also have a question on the output of the MT. The output is supposed > to be signed 32-bit integers, there should be on problem in using > negative integers as seeds for the initialization subroutine right? > > Here's the code snip > Any advice would be sincerely appreciated. Keep up the good work! > > Manuel I'd never seen Mersenne's twister in fortran before as opposed to C. (That would make it naTive around here.) It does eventually compile in fixed format after the line fold-overs are taken care of. Apparently, one can call with an integer to subroutine init_genrand(s) or an array to subroutine init_by_array(init_key,key_length) .. The author warns against 4 functions that are non-standard that have become standard in the meantime: c This program uses the following non-standard intrinsics. c ishft(i,n): If n>0, shifts bits in i by n positions to left. c If n<0, shifts bits in i by n positions to right. c iand (i,j): Performs logical AND on corresponding bits of i and j c ior (i,j): Performs inclusive OR on corresponding bits of i and j c ieor (i,j): Performs exclusive OR on corresponding bits of i and j Given that I have only one node, does it make any sense for me to pursue this? -- "Life in Lubbock, Texas, taught me two things: One is that God loves you and you're going to burn in hell. The other is that sex is the most awful, filthy thing on earth and you should save it for someone you love." ~~ Butch Hancock |