Project CMSSW displayed by LXR |
|
|
|
|||
File indexing completed on 2024-04-06 12:13:24
0001 CCC 0002 CCC from http://cernlib.web.cern.ch/cernlib/download/2005_source/src/mathlib/gen/v/rm48.F 0003 CCC 0004 * 0005 * 0006 * Revision 1.2 2008/04/09 14:11:58 marafino 0007 * Insert new and changed files to accomodate running this code as a producer 0008 * or as a source. 0009 * 0010 * Revision 1.1 2007/02/06 14:40:25 eperez 0011 * first version 0012 * 0013 * Revision 1.2 1996/12/12 16:32:06 cernlib 0014 * Variables ONE and ZERO added to SAVE statement, courtesy R.Veenhof 0015 * 0016 * Revision 1.1.1.1 1996/04/01 15:02:55 mclareni 0017 * Mathlib gen 0018 * 0019 * 0020 SUBROUTINE RM48(RVEC,LENV) 0021 C Double-precision version of 0022 C Universal random number generator proposed by Marsaglia and Zaman 0023 C in report FSU-SCRI-87-50 0024 C based on RANMAR, modified by F. James, to generate vectors 0025 C of pseudorandom numbers RVEC of length LENV, where the numbers 0026 C in RVEC are numbers with at least 48-bit mantissas. 0027 C Input and output entry points: RM48IN, RM48UT. 0028 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 0029 C!!! Calling sequences for RM48: ++ 0030 C!!! CALL RM48 (RVEC, LEN) returns a vector RVEC of LEN ++ 0031 C!!! 64-bit random floating point numbers between ++ 0032 C!!! zero and one. ++ 0033 C!!! CALL RM48IN(I1,N1,N2) initializes the generator from one ++ 0034 C!!! 64-bit integer I1, and number counts N1,N2 ++ 0035 C!!! (for initializing, set N1=N2=0, but to restart ++ 0036 C!!! a previously generated sequence, use values ++ 0037 C!!! output by RM48UT) ++ 0038 C!!! CALL RM48UT(I1,N1,N2) outputs the value of the original ++ 0039 C!!! seed and the two number counts, to be used ++ 0040 C!!! for restarting by initializing to I1 and ++ 0041 C!!! skipping N2*100000000+N1 numbers. ++ 0042 C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 0043 C for 32-bit machines, use IMPLICIT DOUBLE PRECISION 0044 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 0045 DIMENSION RVEC(*) 0046 COMMON/R48ST1/U(97),C,I97,J97 0047 PARAMETER (MODCNS=1000000000) 0048 SAVE CD, CM, TWOM24, NTOT, NTOT2, IJKL,TWOM49, ONE, ZERO 0049 DATA NTOT,NTOT2,IJKL/-1,0,0/ 0050 C 0051 IF (NTOT .GE. 0) GO TO 50 0052 C 0053 C Default initialization. User has called RM48 without RM48IN. 0054 IJKL = 54217137 0055 NTOT = 0 0056 NTOT2 = 0 0057 KALLED = 0 0058 GO TO 1 0059 C 0060 ENTRY RM48IN(IJKLIN, NTOTIN,NTOT2N) 0061 C Initializing routine for RM48, may be called before 0062 C generating pseudorandom numbers with RM48. The input 0063 C values should be in the ranges: 0<=IJKLIN<=900 OOO OOO 0064 C 0<=NTOTIN<=999 999 999 0065 C 0<=NTOT2N<<999 999 999! 0066 C To get the standard values in Marsaglia's paper, IJKLIN=54217137 0067 C NTOTIN,NTOT2N=0 0068 IJKL = IJKLIN 0069 NTOT = MAX(NTOTIN,0) 0070 NTOT2= MAX(NTOT2N,0) 0071 KALLED = 1 0072 C always come here to initialize 0073 1 CONTINUE 0074 IJ = IJKL/30082 0075 KL = IJKL - 30082*IJ 0076 I = MOD(IJ/177, 177) + 2 0077 J = MOD(IJ, 177) + 2 0078 K = MOD(KL/169, 178) + 1 0079 L = MOD(KL, 169) 0080 WRITE(6,'(A,I10,2X,2I10)') ' RM48 INITIALIZED:',IJKL,NTOT,NTOT2 0081 CCC PRINT '(A,4I10)', ' I,J,K,L= ',I,J,K,L 0082 ONE = 1. 0083 HALF = 0.5 0084 ZERO = 0. 0085 DO 2 II= 1, 97 0086 S = 0. 0087 T = HALF 0088 DO 3 JJ= 1, 48 0089 M = MOD(MOD(I*J,179)*K, 179) 0090 I = J 0091 J = K 0092 K = M 0093 L = MOD(53*L+1, 169) 0094 IF (MOD(L*M,64) .GE. 32) S = S+T 0095 3 T = HALF*T 0096 2 U(II) = S 0097 TWOM49 = T 0098 TWOM24 = ONE 0099 DO 4 I24= 1, 24 0100 4 TWOM24 = HALF*TWOM24 0101 C = 362436.*TWOM24 0102 CD = 7654321.*TWOM24 0103 CM = 16777213.*TWOM24 0104 I97 = 97 0105 J97 = 33 0106 C Complete initialization by skipping 0107 C (NTOT2*MODCNS + NTOT) random numbers 0108 DO 45 LOOP2= 1, NTOT2+1 0109 NOW = MODCNS 0110 IF (LOOP2 .EQ. NTOT2+1) NOW=NTOT 0111 IF (NOW .GT. 0) THEN 0112 WRITE(6,'(A,I15)') ' RM48IN SKIPPING OVER ',NOW 0113 DO 40 IDUM = 1, NTOT 0114 UNI = U(I97)-U(J97) 0115 IF (UNI .LT. ZERO) UNI=UNI+ONE 0116 U(I97) = UNI 0117 I97 = I97-1 0118 IF (I97 .EQ. 0) I97=97 0119 J97 = J97-1 0120 IF (J97 .EQ. 0) J97=97 0121 C = C - CD 0122 IF (C .LT. ZERO) C=C+CM 0123 40 CONTINUE 0124 ENDIF 0125 45 CONTINUE 0126 IF (KALLED .EQ. 1) RETURN 0127 C 0128 C Normal entry to generate LENV random numbers 0129 C Modified 13 Marcch '08 to use the CLHEP engines 0130 50 CONTINUE 0131 DO 100 IVEC= 1, LENV 0132 * UNI = U(I97)-U(J97) 0133 * IF (UNI .LT. ZERO) UNI=UNI+ONE 0134 * U(I97) = UNI 0135 * I97 = I97-1 0136 * IF (I97 .EQ. 0) I97=97 0137 * J97 = J97-1 0138 * IF (J97 .EQ. 0) J97=97 0139 * C = C - CD 0140 * IF (C .LT. ZERO) C=C+CM 0141 * UNI = UNI-C 0142 * IF (UNI .LT. ZERO) UNI=UNI+ONE 0143 * RVEC(IVEC) = UNI 0144 UNI = BHGPYR(IDUMMY) 0145 RVEC(IVEC) = UNI 0146 C Replace exact zeros by 2**-49 0147 * IF (UNI .EQ. ZERO) THEN 0148 * RVEC(IVEC) = TWOM49 0149 * ENDIF 0150 100 CONTINUE 0151 NTOT = NTOT + LENV 0152 IF (NTOT .GE. MODCNS) THEN 0153 NTOT2 = NTOT2 + 1 0154 NTOT = NTOT - MODCNS 0155 ENDIF 0156 RETURN 0157 C Entry to output current status 0158 ENTRY RM48UT(IJKLUT,NTOTUT,NTOT2T) 0159 IJKLUT = IJKL 0160 NTOTUT = NTOT 0161 NTOT2T = NTOT2 0162 RETURN 0163 END
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.2.1 LXR engine. The LXR team |
|