1 SUBROUTINE initgn(isdtyp) 2C********************************************************************** 3C 4C SUBROUTINE INITGN(ISDTYP) 5C INIT-ialize current G-e-N-erator 6C 7C Reinitializes the state of the current generator 8C 9C This is a transcription from Pascal to Fortran of routine 10C Init_Generator from the paper 11C 12C L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package 13C with Splitting Facilities." ACM Transactions on Mathematical 14C Software, 17:98-111 (1991) 15C 16C 17C Arguments 18C 19C 20C ISDTYP -> The state to which the generator is to be set 21C 22C ISDTYP = -1 => sets the seeds to their initial value 23C ISDTYP = 0 => sets the seeds to the first value of 24C the current block 25C ISDTYP = 1 => sets the seeds to the first value of 26C the next block 27C 28C INTEGER ISDTYP 29C 30C********************************************************************** 31C .. Parameters .. 32 INTEGER*4 numg 33 PARAMETER (numg=32) 34C .. 35C .. Scalar Arguments .. 36 INTEGER*4 isdtyp 37C .. 38C .. Scalars in Common .. 39 INTEGER*4 a1,a1vw,a1w,a2,a2vw,a2w,m1,m2 40C .. 41C .. Arrays in Common .. 42 INTEGER*4 cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg), 43 + lg2(numg) 44 LOGICAL qanti(numg) 45C .. 46C .. Local Scalars .. 47 INTEGER*4 g 48C .. 49C .. External Functions .. 50 LOGICAL qrgnin 51 INTEGER*4 mltmod 52 EXTERNAL qrgnin,mltmod 53C .. 54C .. External Subroutines .. 55 EXTERNAL getcgn 56C .. 57C .. Common blocks .. 58 COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1, 59 + cg2,qanti 60C .. 61C .. Save statement .. 62 SAVE /globe/ 63C .. 64C .. Executable Statements .. 65C Abort unless random number generator initialized 66 IF (qrgnin()) GO TO 10 67 WRITE (*,*) ' INITGN called before random number generator ', 68 + ' initialized -- abort!' 69 CALL XSTOPX 70 + (' INITGN called before random number generator initialized') 71 72 10 CALL getcgn(g) 73 IF ((-1).NE. (isdtyp)) GO TO 20 74 lg1(g) = ig1(g) 75 lg2(g) = ig2(g) 76 GO TO 50 77 78 20 IF ((0).NE. (isdtyp)) GO TO 30 79 CONTINUE 80 GO TO 50 81C do nothing 82 30 IF ((1).NE. (isdtyp)) GO TO 40 83 lg1(g) = mltmod(a1w,lg1(g),m1) 84 lg2(g) = mltmod(a2w,lg2(g),m2) 85 GO TO 50 86 87 40 CALL XSTOPX ('ISDTYP NOT IN RANGE') 88 89 50 cg1(g) = lg1(g) 90 cg2(g) = lg2(g) 91 RETURN 92 93 END 94