00001 SUBROUTINE RANMAR(RVEC,LENV)
00002 C ----------------------------------------------------------------------
00003 C<<<<<FUNCTION RANMAR(IDUMM)
00004 C CERNLIB V113, VERSION WITH AUTOMATIC DEFAULT INITIALIZATION
00005 C Transformed to SUBROUTINE to be as in CERNLIB
00006 C AM.Lutz November 1988, Feb. 1989
00007 C
00008 C!Universal random number generator proposed by Marsaglia and Zaman
00009 C in report FSU-SCRI-87-50
00010 C modified by F. James, 1988 and 1989, to generate a vector
00011 C of pseudorandom numbers RVEC of length LENV, and to put in
00012 C the COMMON block everything needed to specify currrent state,
00013 C and to add input and output entry points RMARIN, RMARUT.
00014 C
00015 C Unique random number used in the program
00016 C ----------------------------------------------------------------------
00017 COMMON / INOUT / INUT,IOUT
00018 DIMENSION RVEC(*)
00019 COMMON/RASET1/U(97),C,I97,J97
00020 PARAMETER (MODCNS=1000000000)
00021 DATA NTOT,NTOT2,IJKL/-1,0,0/
00022 C
00023 IF (NTOT .GE. 0) GO TO 50
00024 C
00025 C Default initialization. User has called RANMAR without RMARIN.
00026 IJKL = 54217137
00027 NTOT = 0
00028 NTOT2 = 0
00029 KALLED = 0
00030 GO TO 1
00031 C
00032 ENTRY RMARIN(IJKLIN, NTOTIN,NTOT2N)
00033 C Initializing routine for RANMAR, may be called before
00034 C generating pseudorandom numbers with RANMAR. The input
00035 C values should be in the ranges: 0<=IJKLIN<=900 OOO OOO
00036 C 0<=NTOTIN<=999 999 999
00037 C 0<=NTOT2N<<999 999 999!
00038 C To get the standard values in Marsaglia-s paper, IJKLIN=54217137
00039 C NTOTIN,NTOT2N=0
00040 IJKL = IJKLIN
00041 NTOT = MAX(NTOTIN,0)
00042 NTOT2= MAX(NTOT2N,0)
00043 KALLED = 1
00044 C always come here to initialize
00045 1 CONTINUE
00046 IJ = IJKL/30082
00047 KL = IJKL - 30082*IJ
00048 I = MOD(IJ/177, 177) + 2
00049 J = MOD(IJ, 177) + 2
00050 K = MOD(KL/169, 178) + 1
00051 L = MOD(KL, 169)
00052 WRITE(IOUT,201) IJKL,NTOT,NTOT2
00053 201 FORMAT(1X,' RANMAR INITIALIZED: ',I10,2X,2I10)
00054 DO 2 II= 1, 97
00055 S = 0.
00056 T = .5
00057 DO 3 JJ= 1, 24
00058 M = MOD(MOD(I*J,179)*K, 179)
00059 I = J
00060 J = K
00061 K = M
00062 L = MOD(53*L+1, 169)
00063 IF (MOD(L*M,64) .GE. 32) S = S+T
00064 3 T = 0.5*T
00065 2 U(II) = S
00066 TWOM24 = 1.0
00067 DO 4 I24= 1, 24
00068 4 TWOM24 = 0.5*TWOM24
00069 C = 362436.*TWOM24
00070 CD = 7654321.*TWOM24
00071 CM = 16777213.*TWOM24
00072 I97 = 97
00073 J97 = 33
00074 C Complete initialization by skipping
00075 C (NTOT2*MODCNS + NTOT) random numbers
00076 DO 45 LOOP2= 1, NTOT2+1
00077 NOW = MODCNS
00078 IF (LOOP2 .EQ. NTOT2+1) NOW=NTOT
00079 IF (NOW .GT. 0) THEN
00080 WRITE (IOUT,'(A,I15)') ' RMARIN SKIPPING OVER ',NOW
00081 DO 40 IDUM = 1, NTOT
00082 UNI = U(I97)-U(J97)
00083 IF (UNI .LT. 0.) UNI=UNI+1.
00084 U(I97) = UNI
00085 I97 = I97-1
00086 IF (I97 .EQ. 0) I97=97
00087 J97 = J97-1
00088 IF (J97 .EQ. 0) J97=97
00089 C = C - CD
00090 IF (C .LT. 0.) C=C+CM
00091 40 CONTINUE
00092 ENDIF
00093 45 CONTINUE
00094 IF (KALLED .EQ. 1) RETURN
00095 C
00096 C Normal entry to generate LENV random numbers
00097 50 CONTINUE
00098 DO 100 IVEC= 1, LENV
00099 UNI = U(I97)-U(J97)
00100 IF (UNI .LT. 0.) UNI=UNI+1.
00101 U(I97) = UNI
00102 I97 = I97-1
00103 IF (I97 .EQ. 0) I97=97
00104 J97 = J97-1
00105 IF (J97 .EQ. 0) J97=97
00106 C = C - CD
00107 IF (C .LT. 0.) C=C+CM
00108 UNI = UNI-C
00109 IF (UNI .LT. 0.) UNI=UNI+1.
00110 C Replace exact zeroes by uniform distr. *2**-24
00111 IF (UNI .EQ. 0.) THEN
00112 UNI = TWOM24*U(2)
00113 C An exact zero here is very unlikely, but lets be safe.
00114 IF (UNI .EQ. 0.) UNI= TWOM24*TWOM24
00115 ENDIF
00116 RVEC(IVEC) = UNI
00117 100 CONTINUE
00118 NTOT = NTOT + LENV
00119 IF (NTOT .GE. MODCNS) THEN
00120 NTOT2 = NTOT2 + 1
00121 NTOT = NTOT - MODCNS
00122 ENDIF
00123 RETURN
00124 C Entry to output current status
00125 ENTRY RMARUT(IJKLUT,NTOTUT,NTOT2T)
00126 IJKLUT = IJKL
00127 NTOTUT = NTOT
00128 NTOT2T = NTOT2
00129 RETURN
00130 END
