c Code module file pranf.F c c----------------------------------------------------------------------- c c This module contains the parallel pseudorandom number generator Pranf. c It is designed to duplicate the Cray PRNG "ranf" in a portable and c wordlength independent form, compatible with basic Fortran-77, and to c generate multiple streams of random 48-bit floating point values c between zero and one. Each stream is autonomous, independent, and c requires its own seed. 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, StartVal, Seeds ) c 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, StartVal, 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 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 c----------------------------------------------------------------------- c======================================================================= c c P A R A L L E L 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 ) 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 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, StartVal, 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 StartVal 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, StartVal is set to a special c value, [3281, 4041, 595, 2376], 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, StartVal, 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, StartVal 10 format( ' Rans: NIn = ', i10, ' StartVal = ', 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( StartVal .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( StartVal ) 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=======================================================================