On 2008-04-24 11:45:55 -0300, mjm2114@[EMAIL PROTECTED]
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.hiro****ma-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 Ni****mura 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 Ni****mura,
> 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 do***entation 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.hiro****ma-u.ac.jp/~m-mat/MT/emt.html
> c email: m-mat @[EMAIL PROTECTED]
math.sci.hiro****ma-u.ac.jp (remove space)
> c
> c-----------------------------------------------------------------------
> c FORTRAN77 translation by Tsuyo**** 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, ****fts bits in i by n positions to left.
> c If n<0, ****fts 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
|
