Prev: Can I create header file dependent by gfortran
Next: Bug in gfortran when reading end of file?
From: mjm2114 on 24 Apr 2008 10:45
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

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: mjm2114 on 24 Apr 2008 11:22
Thanks for your reply. Maybe I wasn't clear. I take the first block of
624 outputs from the master node to seed the first worker node. I then
take the second block of 624 outputs from the master node to seed the
second worker node, and so forth. In this way, every worker node has
been initialized with a different array of seeds. The master node is
only used to produce seeds for the worker nodes, it is not used for
anything else.
From: mjm2114 on 24 Apr 2008 12:39
On Apr 24, 12:01 pm, Gordon Sande <g.sa...(a)worldnet.att.net> wrote:
> On 2008-04-24 12:22:49 -0300, mjm2...(a)columbia.edu said:
>
> > Thanks for your reply. Maybe I wasn't clear. I take the first block of
> > 624 outputs from the master node to seed the first worker node. I then
> > take the second block of 624 outputs from the master node to seed the
> > second worker node, and so forth. In this way, every worker node has
> > been initialized with a different array of seeds. The master node is
> > only used to produce seeds for the worker nodes, it is not used for
> > anything else.
>
> Your description still makes it sound like each node will be only
> a delayed version of a lower numbered node.
>
> The internal state of the PRNG is usually called the seed. There may be
> shortcut initial startups, labelled seeds, for those generators that have
> large seeds like MT. Leads to confusion on the meaning of the words when
> the adjectives are dropped. The purpose of a full seed is to pick up again
> EXACTLY where you left off. Having the master step along and give seeds
> to each slave just means that the slaves will exactly follow what the
> master would have done. Exact reproducability of PRNGs is important for
> debugging. That is what the saving and restoring of seeds is all about.
>
> It would be easy to form the conclusion that you are rather confused
> about what the role of the seed is. Neither your initial nor followup
> posting serve to dispell that. For MT, "THE SEED" is all of the 600+
> words. It is not an array of seeds. How else do you think they achieve
> their immense cycle length?

Thanks again for your comments. You are correct, my terminology was
confusing.
If we consider "seed" as synonymous with internal state then MT has a
seed of 624 32-bit integers.
To initialize the internal state of an MT in each worker node I use a
contiguous chunk of outputs from the MT in the master node. This chunk
happens to have length 624 (recommended by MT literature) but is is
not going to be the same as the initial internal state in the worker
MT (algorithm takes care of this), nor is it equal to the previous
internal state in the master node (because this previous state was
only used to create the last number in the chunk). Hope this clears
things up.

thanks again for your comments

Manuel
From: e p chandler on 24 Apr 2008 20:25
On Apr 24, 3:52 pm, Gordon Sande <g.sa...(a)worldnet.att.net> wrote:

> In the case of the "one word internal state" a case could be made
> that the seed could take states that were not internally possible
> but that the reported internal state would be a seed as well. This
> was for the situation where the 31 bits had to end in "01" (3 mod 8)
> so some of bits in a 32 bit seed could be invalid. Seed is a curious
> uasge at best but seems to be rather well establshed.

True only if the additive constant is ZERO.

> The OP seemed to be saying that each word in the many word internal
> state was a seed so in an attempt to draw his attention to the whole
> thing I used "full seed" on the fly. Better suggestions are welcome. ;-)

I'll stay away from this. :-(. IMO it's more important to make sure
that ALL bits of the initial state are initialized AND that all bits
are initialized properly. [For a counter example look at how MS
Basic's RANDOMIZE n works in version 5.x]

[Sorry it's been years since I've read about GFSR routines. Back then
not all the bits could be zero. Now it seems that an excess of zero
bits in the intial state is bad.... Too complex for me to understand
anyway. :-(.]

-- e

 | 
Pages: 1
Prev: Can I create header file dependent by gfortran
Next: Bug in gfortran when reading end of file?