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