c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 program RanfTest c======================================================================= c----------------------------------------------------------------------- c c This program does test runs of the parallel pseudorandom number gener- c ator Pranf. It is designed to generate multiple streams of random 48- c bit floating point values between zero and one. Each stream is auto- c nomous, independent, and requires its own seed. it is based on the c Cray RNG known as "ranf". The method for generating multiple streams c by partitioning the ranf sequence were provided by Dr. Mark J. Durst, c of Lawrence Livermore National Laboratory. c c The original version of Pranf was written in SISAL for the parallel c simulated annealing program PSA, written in SISAL. This fortran ver- c sion was written as part of a comparative performance study, by T. M. c DeBoni in conjunction with John Feo and David Cann, all of Lawrence c Livermore National Laboratory, c c Here's how to use these functions: c c Initialization: call Rans( N, Seed1, Seeds ) c c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c This returns an array of M seeds, where M = N if N is odd, or N+1 c otherwise; the returned seeds have values that partition the basic c ranf cycle of 2**46 pseudorandom reals in [0,1] into independent sub- c sequences. The second argument, Seed1, can be a cycle-starting seed c for the full-period generator; if it is zero, a special seed will be c used that produces statistically desirable behavior. c c Use: call Ranf( Seed, RandNum ) c c This returns a pseudorandom real and a seed to be passed back in the c next invocation. The returned seed carries the state normally hidden c in imperative generators. c c Note: the Seed returned by Ranf is actually the next RandNum, c encoded in a special form, to be described below. c c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c The individual subsequences each carry the statistical properties of c the c whole cycle, and this is a good generator, adapted from that c used at LLNL on the Cray computers by the fortran compilers. However, c care must be taken that the subsequences are used consistently. If c processes trade seeds, undesirable statistical behavior may result. c c Note that a special data structure is needed for this generator; each c seed returned by the generator is really an array of four integers. c These are intended to be the four digits of a radix-2**12 number. All c calculations internal to the generator and partitioner are carried c out in piecewise fashion, using these numbers. This allows the gener- c ator to be easily duplicated on machines of varying word length. The c random number returned is calculated as if a 48-bit mantissa were c the desired representation. This is consistent with Cray numerics. On c machines with shorter words, the number's lowest bits will be lost. c However, it should be possible to reprogram this generator to produce c a different final representation for the random number, if desired. c c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c----------------------------------------------------------------------- c======================================================================= c c C O N S T A N T D E F I N I T I O N S c c======================================================================= c c Following are sample contents of the file "pran.parm", which is c included in each routine. It contains parameter definitions that are c used to size the data structures. The include statement is of the form #include "pranf.parm" c c c**** Program control and storage sizing parameters: c**** ---------------------------------------------- c parameter c**** The upper limits possible in this compilation: c 1 ( MaxParallelism = 128, c**** Some array sizing values: c 2 MaxRandNumStreams = MaxParallelism, c 3 IRandNumSize = 4, IBinarySize = 48, c 4 Mod4096DigitSize = 12, c 5 IAmountofParallelism = 128 ) c**** Data structuring parameters: c**** --------------------------- c**** Other control values: c**** -------------------- c parameter ( IZero = 0, Inunit = 5, IOutunit = 6, c 1 IRstunit = 7, MyTrue = 1, MyFalse = 0 ) c parameter ( MaxNumberofNumbers = 100000 ) c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c======================================================================= c c S T O R A G E D E F I N I T I O N S c c======================================================================= c c First, some discussion of the data sturctures we use in this program. c Since Fortran has no data structuring abilities except arrays, all of c these structures are represented by arrays. Special access and modifi- c cation c routines are provied to maintain the appropriate data abs- c tractions. These routines are given elsewhere. c c**** This is one seed used in the random number generator; it is repre- c**** sented as an array of four modulo-4096 (12-bit) integers: c Seed( IRandNumSize ) c c**** This is the array of seeds used in the parallel random number c**** generator: c Seeds( MaxRandNumStreams, IRandNumSize ) c c**** Here, now, are the actual storage elements used in the main program. c**** Others of those mentioned above will be found in other routines. c integer NumberofStreams, SeedIn, NumberofRandNums, Seed, Seeds c**** Added array below on 05 Mar 98 - TMD integer AllSeeds real*8 RandNum, AC, ACSQ, RandNums dimension 1 RandNums( MaxNumberofNumbers ), c**** Added array below on 05 Mar 98 - TMD * AllSeeds( MaxNumberofNumbers, IRandNumSize ), 2 AC( MaxRandNumStreams ), 3 ACSQ( MaxRandNumStreams ), 4 Seeds( MaxRandNumStreams, IRandNumSize ), 5 Seed( IRandNumSize ) c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c======================================================================= c c M A I N P R O G R A M c c======================================================================= c**** Set up input and output files. call FileSetup c**** Read the input data. call GetInput( NumberofStreams, SeedIn, NumberofRandNums ) if ( NumberofStreams .gt. MaxRandNumStreams ) then write( *, 60 ) NumberofStreams 60 format( ' NumberofStreams = ', i10, ' > MaxRandNumStreams; ' ) endif if ( IAmountofParallelism .gt. MaxParallelism ) then write( *, 70 ) IAmountofParallelism 70 format( ' IAmountofParallelism = ', i10, ' > MaxParallelism; ' ) endif if ( NumberofRandNums .gt. MaxNumberofNumbers ) then write( *, 80 ) 80 format( ' NumberofRandNums = ', i10, ' > MaxNumberofNumbers; ' ) endif if( ( NumberofStreams .gt. MaxRandNumStreams ) .or. 1 ( IAmountofParallelism .gt. MaxParallelism ) .or. 2 ( NumberofRandNums .gt. MaxNumberofNumbers ) ) then write( *, 90 ) 90 format( ' exiting. ' ) stop endif c**** Initialize the parallel random number generator to the number of c**** separate streams specified by the amount of parallelism needed. call Rans( NumberofStreams, SeedIn, Seeds ) c**** Generate a stream of numbers. do 110 j = 1, NumberofStreams AC(j) = 0d0 ACSQ(j) = 0d0 call GetSeed( Seeds, j, Seed ) do 100 I = 1, NumberofRandNums call Ranf( Seed, RandNum ) RandNums( I ) = RandNum AC(j) = BSJ0(RandNum) ACSQ(j) = BSK0(RandNum) c**** Added data structure below on 05 Mar 98 - TMD AllSeeds( I, 1 ) = Seed(1) AllSeeds( I, 2 ) = Seed(2) AllSeeds( I, 3 ) = Seed(3) AllSeeds( I, 4 ) = Seed(4) AC(j) = AC(j) + RandNum ACSQ(j) = ACSQ(j) + RandNum * RandNum 100 continue 110 continue c**** Write out the data. write( *, 125 ) NumberofRandNums 125 format( ' Stream of ', i10, ' random numbers:' ) c**** changed output on 05 Mar 98 - TMD write( *, 150 ) ( ((AllSeeds(I, J), J = 1, 4), 1 (RandNums( I ), RandNums(I))), 2 I = 1, NumberofRandNums ) c150 format( 3d22.15 ) 150 format( 4(i4, x), z16, x, g22.15 ) write( *, 155 ) 155 format( ' Stats: ' ) do 170 j = 1, NumberofStreams write( *, 160 ) j, AC(j), j, ACSQ(j) 160 format( ' AC(',i3,') = ', d22.15 ' ACSQ(',i3,' = ', d22.15 ) 170 continue c**** Close all the files. call CloseAllFiles stop end c end Main c======================================================================= c c D A T A S T R U C T U R E M A N I P U L A T I O N C O D E c c======================================================================= subroutine PutSeed( Seeds, Seed, SeedIndex ) c----------------------------------------------------------------------- #include "pranf.parm" integer Seeds, Seed, SeedIndex dimension Seeds( MaxRandNumStreams, IRandNumSize ), 2 Seed( IRandNumSize ) c----------------------------------------------------------------------- c write( *, 10 ) 10 format( ' PutSeed', /, '-------' ) Seeds( SeedIndex, 1 ) = Seed( 1 ) Seeds( SeedIndex, 2 ) = Seed( 2 ) Seeds( SeedIndex, 3 ) = Seed( 3 ) Seeds( SeedIndex, 4 ) = Seed( 4 ) c write( *, 20 ) 20 format( ' leaving PutSeed' ) return end c end PutSeed c======================================================================= subroutine GetSeed( Seeds, SeedIndex, Seed ) c----------------------------------------------------------------------- #include "pranf.parm" integer Seeds, SeedIndex, Seed dimension Seeds( MaxRandNumStreams, IRandNumSize ), 2 Seed( IRandNumSize ) c----------------------------------------------------------------------- c write( *, 10 ) 10 format( ' GetSeed', /, '-------' ) Seed( 1 ) = Seeds( SeedIndex, 1 ) Seed( 2 ) = Seeds( SeedIndex, 2 ) Seed( 3 ) = Seeds( SeedIndex, 3 ) Seed( 4 ) = Seeds( SeedIndex, 4 ) c write( *, 20 ) 20 format( ' leaving GetSeed' ) return end c end GetSeed c======================================================================= c c A U X I L I A R Y R O U T I N E S c c======================================================================= subroutine CloseAllFiles c----------------------------------------------------------------------- #include "pranf.parm" c----------------------------------------------------------------------- c write( *, 10 ) 10 format( ' CloseAllFiles ', /, '--------------' ) c close( IOutunit ) return end c end CloseAllFiles c======================================================================= subroutine FileSetup c----------------------------------------------------------------------- #include "pranf.parm" character*23 InputFileName, OutputFileName c----------------------------------------------------------------------- c write( *, 10 ) c10 format( ' Enter name of file containing data tuples: ' ) c read( *, 20 ) InputFileName c20 format( A ) c write( *, 30 ) c30 format( ' Enter name of file to receive all text output: ' ) c read( *, 40 ) OutputFileName c40 format( A ) c open( file=InputFileName, unit=Inunit ) c open( file=OutputFileName, unit=IOutunit ) c write( *, 50 ) 50 format( ' FileSetup', /, '----------' ) return end c end FileSetup c======================================================================= subroutine GetInput( NumberofStreams, SeedIn, NumberofRandNums ) c----------------------------------------------------------------------- #include "pranf.parm" integer NumberofStreams, SeedIn, NumberofRandNums c----------------------------------------------------------------------- write( *, 10 ) 10 format( ' GetInput', /, '---------' ) c**** Read the control values. read( *, 20 ) NumberofStreams, SeedIn, NumberofRandNums 20 format( 3i10 ) c**** Echo the control values. write( *, 110 ) NumberofStreams, SeedIn, NumberofRandNums 110 format( ' INPUT DATA:', /, '===========', /, 1 ' NumberofStreams = ', i10, ' SeedIn = ', i10, 2 ' NumberofRandNums = ', i10, / ) write( *, 210 ) 210 format( 'Closing input file.', / ) c**** All done with the input file. c close( Inunit ) return end c end GetInput c======================================================================= subroutine PutOutput( ) c----------------------------------------------------------------------- #include "pranf.parm" c----------------------------------------------------------------------- write( *, 100 ) 100 format( ' OUTPUT DATA:', /, '============' ) return end c end PutOutput c======================================================================= c c R A N D O M N U M B E R G E N E R A T O R c c======================================================================= subroutine Ranf( Seed, RandNum ) c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c----------------------------------------------------------------------- c c This is the one of the two user-callable routines of a linear congru- c ential random number generator, modeled after the RANF generator for c the Cray X/MP. This routine generates the next random number in the c sequence and a new seed for the remainder of the sequence.The seed and c the random number are the same, but are returned in different form: c the random number is a fortran 'real', but the seed is an array of c four integers, each containing an 12-bit value that is used internal- c ly to the generator as one of four digits of a radix-4096 number. c c It returns the new random number and a new seed. c c----------------------------------------------------------------------- #include "pranf.parm" integer Seed real*8 RandNum c*********************************************************************** c**** Data common to the PRanf package. integer Multiplier, DefaultSeed real*8 Divisor dimension Multiplier( IRandNumSize ), DefaultSeed( IRandNumSize ), 1 Divisor( IRandNumSize ) dimension Seed( IRandNumSize ) common / PRanf / Multiplier, DefaultSeed, Divisor data Multiplier / 373, 3707, 1442, 647 /, 1 DefaultSeed / 3281, 4041, 595, 2376 / , c**** changing default seed as per MJD 05 Mar 98 - TMD c**** changing default seed back to original one, above. c 1 DefaultSeed / 1645, 4065, 1676, 700 / , 2 Divisor / 281474976710656.0, 68719476736.0, 16777216.0, 3 4096.0 / c**** End of PRanf common data c*********************************************************************** c----------------------------------------------------------------------- c write( *, 10 ) 10 format( ' Ranf', /, '-----' ) RandNum = float( Seed( 4 ) ) / Divisor( 4 ) + 1 float( Seed( 3 ) ) / Divisor( 3 ) + 2 float( Seed( 2 ) ) / Divisor( 2 ) + 3 float( Seed( 1 ) ) / Divisor( 1 ) call ranfmodmult( Multiplier, Seed, Seed ) c write( *, 20 ) 20 format( ' leaving Ranf' ) return end c end Ranf c======================================================================= subroutine Rans( NIn, Seed1, SeedArray ) c----------------------------------------------------------------------- c c This routine divides the sequence of random numbers into N >= NIn sub- c sequences, each with its own seed. The seeds for the independent sub- c sequences are returned in eedArray. If Seed1 is zero, all zeroes will c be returned as a result of the modulus arithmetic involved in parti- c tioning the main sequence. To prevent this, Seed1 is set to a special c value, [3281, 4041, 595, 2376], which is statistically a good starting c changed line above to line below - 05 mar 98 - TMD c value, [1645, 4065, 1676, 700], which is statistically a good starting c seed. The main sequence is then divided into the N subsequences (where c N is odd and >= NIn) by dividing its period (2**46) by N. A prime c value for N is best, but an odd value is better than an even one. c c Then, seed(i) = seed(i-1) * (a**k mod 2**48), and 1<=k<=N, c c where, 'a' is [373, 3707, 1442, 647], the multiplier used by the LCG c which is being partitioned. c c The number of streams must be odd; if NIn is even N will be NIn+1, c and an extra stream of random numbers will be available that may not c be c used. c c It returns an array of seeds, each an array of 4 integers that are c used as c the digits of a four-digit modulo-4096 integer. c c----------------------------------------------------------------------- #include "pranf.parm" integer NIn, Seed1, DefaultValues, SeedArray, N, atothek, K, 1 KBinary, InSeed, OutSeed c*********************************************************************** c**** Data common to the PRanf package. integer Multiplier, DefaultSeed real*8 Divisor dimension Multiplier( IRandNumSize ), 1 DefaultSeed( IRandNumSize ), 1 Divisor( IRandNumSize ) common / PRanf / Multiplier, DefaultSeed, Divisor c**** End of PRanf common data c*********************************************************************** dimension DefaultValues( IRandNumSize ), 1 SeedArray( MaxRandNumStreams, IRandNumSize ), 2 atothek( IRandNumSize ), 3 K( IRandNumSize ), 4 KBinary( IBinarySize ), 5 InSeed( IRandNumSize ), 6 OutSeed( IRandNumSize ) c----------------------------------------------------------------------- c**** Statement function the check for argument numbers being odd. NIsOdd( m ) = mod( m, 2 ) write( *, 10 ) NIn, Seed1 10 format( ' Rans: NIn = ', i10, ' Seed1 = ', i10, /, 1 '------' ) c**** Make sure we are generating an odd number of random number c**** sequences. if( NisOdd( NIn ) .eq. 1 ) then N = NIN else N = NIn + 1 endif c**** Set up the initial seed to either a legal input value or its c**** default value sets. The input integer, if nonzero, is used for c**** the first of the four modulo-4096 digits of the actual initial c**** seed. c**** The First element of the returned SeedArray will be used here, c**** since at least one seed will be returned, and the first seed of c**** the set returned will be the seed for the entire wheel. if( Seed1 .eq. IZero ) then SeedArray( 1, 1 ) = DefaultSeed( 1 ) SeedArray( 1, 2 ) = DefaultSeed( 2 ) SeedArray( 1, 3 ) = DefaultSeed( 3 ) SeedArray( 1, 4 ) = DefaultSeed( 4 ) else SeedArray( 1, 1 ) = abs( Seed1 ) SeedArray( 1, 2 ) = IZero SeedArray( 1, 3 ) = IZero SeedArray( 1, 4 ) = IZero endif c**** 'a' is the multiplier for the Ranf linear congruential generator. if( N .eq. 1 ) then c******* If only one stream of random numbers is needed, do not bother c******* to raise 'a' to the first power. atothek( 1 ) = Multiplier( 1 ) atothek( 2 ) = Multiplier( 2 ) atothek( 3 ) = Multiplier( 3 ) atothek( 4 ) = Multiplier( 4 ) else c******* more than one stream is needed; generate the Kth seed by c******* multiplying the K-1st seed by the multiplier raised to the c******* Nth power. call ranfk( N, K ) call ranfkbinary( K, KBinary ) call ranfatok( Multiplier, KBinary, atothek ) do 100 I = 2, N InSeed( 1 ) = SeedArray( I-1, 1 ) InSeed( 2 ) = SeedArray( I-1, 2 ) InSeed( 3 ) = SeedArray( I-1, 3 ) InSeed( 4 ) = SeedArray( I-1, 4 ) call ranfmodmult( InSeed, atothek, OutSeed ) SeedArray( I, 1 ) = OutSeed( 1 ) SeedArray( I, 2 ) = OutSeed( 2 ) SeedArray( I, 3 ) = OutSeed( 3 ) SeedArray( I, 4 ) = OutSeed( 4 ) 100 continue endif write( *, 110 ) 110 format( ' leaving Rans' ) do 120 i = 1, N write( *, 105 ) i, (SeedArray(i, j), j = 1, IRandNumSize) 105 format( ' SeedArray(', i3, ') = ', 4(i5,x) ) 120 continue return end c end Rans c======================================================================= subroutine ranfatok( a, Kbinary, atothek ) c----------------------------------------------------------------------- c c This routine calculates a to the Kth power, mod 2**48. K is a binary c number. c c It returns the calculated value as an array of four radix-4096 digits. c c----------------------------------------------------------------------- #include "pranf.parm" integer a, KBinary, atothek, asubi dimension a( IRandNumSize ), 1 KBinary( IBinarySize ), 2 atothek( IRandNumSize ), 3 asubi( IRandNumSize ) c----------------------------------------------------------------------- c**** The following amounts to the first iteration of a 46-loop. asubi( 1 ) = a( 1 ) asubi( 2 ) = a( 2 ) asubi( 3 ) = a( 3 ) asubi( 4 ) = a( 4 ) atothek( 1 ) = 1 atothek( 2 ) = IZero atothek( 3 ) = IZero atothek( 4 ) = IZero do 100 I = 1, 45 if( KBinary( I ) .ne. IZero ) then call ranfmodmult( atothek, asubi, atothek ) endif call ranfmodmult( asubi, asubi, asubi ) 100 continue return end c ranfatok c======================================================================= function iranfeven( N ) c----------------------------------------------------------------------- c c This function checks the parity of the argument integer. c c It returns one if the argument is even and zero otherwise. c c----------------------------------------------------------------------- #include "pranf.parm" integer N c----------------------------------------------------------------------- if( mod( N, 2 ) .eq. 0 ) then iranfeven = 1 else iranfeven = 0 endif return end c end iranfeven c======================================================================= subroutine ranfk( N, K ) c----------------------------------------------------------------------- c c This routine calculates 2**46/N, which should be the period of each c of the subsequences of random numbers that are being created. Both the c numerator and the result of this calculation are represented as an c array of four integers, each of which is one digit of a four-digit c radix-4096 number. The numerator is represented as (1024, 0, 0, 0 ), c using base ten digits. c c It returns the result of the division. c c----------------------------------------------------------------------- #include "pranf.parm" integer N, K, nn, r4, r3, r2, r1, q4, q3, q2, q1 dimension K( IRandNumSize ) c----------------------------------------------------------------------- nn = N + iranfeven( N ) q4 = 1024 / nn r4 = 1024 - (nn * q4) q3 = (r4 * 4096) / nn r3 = (r4 * 4096) - (nn * q3) q2 = (r3 * 4096) / nn r2 = (r3 * 4096) - (nn * q2) q1 = (r2 * 4096) / nn K( 1 ) = q1 K( 2 ) = q2 K( 3 ) = q3 K( 4 ) = q4 return end c end ranfk c======================================================================= subroutine ranfkbinary( K, KBinary ) c----------------------------------------------------------------------- c c This routine calculates the binary expansion of the argument K, which c is a 48-bit integer represented as an array of four 12-bit integers. c c It returns an array of 48 binary values. c c----------------------------------------------------------------------- #include "pranf.parm" integer K, KBinary, X, Bits dimension K( IRandNumSize ), 1 KBinary( IBinarySize ), 2 Bits( Mod4096DigitSize ) c----------------------------------------------------------------------- do 300 I = 1, 4 X = K( I ) / 2 Bits( 1 ) = iranfodd( K( I ) ) do 100 J = 2, Mod4096DigitSize Bits( J ) = iranfodd( X ) X = X / 2 100 continue do 200 J = 1, Mod4096DigitSize KBinary( (I-1)*Mod4096DigitSize + J ) = Bits( J ) 200 continue 300 continue return end c end ranfkbinary c======================================================================= subroutine ranfmodmult( A, B, C ) c----------------------------------------------------------------------- c c Ths routine calculates the product of the first two arguments. All c three arguements are represented as arrays of 12-bit integers, each c making up the digits of one radix-4096 number. The multiplication is c done piecemeal. c c It returns the product in the third argument. c c----------------------------------------------------------------------- #include "pranf.parm" integer A, B, C, j1, j2, j3, j4, k1, k2, k3, k4 dimension A( IRandNumSize ), 1 B( IRandNumSize ), 2 C( IRandNumSize ) c----------------------------------------------------------------------- c write( *, 10 ) 10 format( ' ranfmodmult', /, '------------' ) j1 = A( 1 ) * B( 1 ) j2 = A( 1 ) * B( 2 ) + A( 2 ) * B( 1 ) j3 = A( 1 ) * B( 3 ) + A( 2 ) * B( 2 ) + A( 3 ) * B( 1 ) j4 = A( 1 ) * B( 4 ) + A( 2 ) * B( 3 ) + A( 3 ) * B( 2 ) + 1 A( 4 ) * B( 1 ) k1 = j1 k2 = j2 + k1 / 4096 k3 = j3 + k2 / 4096 k4 = j4 + k3 / 4096 C( 1 ) = mod( k1, 4096 ) C( 2 ) = mod( k2, 4096 ) C( 3 ) = mod( k3, 4096 ) C( 4 ) = mod( k4, 4096 ) c write( *, 20 ) 20 format( ' leaving ranfmodmult' ) return end c end ranfmodmult c======================================================================= function iranfodd( N ) c----------------------------------------------------------------------- c c This function checks the parity of the argument integer. c c It returns one if the argument is odd and zero otherwise. c c----------------------------------------------------------------------- #include "pranf.parm" integer N c----------------------------------------------------------------------- if( mod( N, 2 ) .eq. 0 ) then iranfodd = 0 else iranfodd = 1 endif return end c end iranfodd c=======================================================================