Mercurial > octave
changeset 3188:df7c57a6639d
[project @ 1998-10-15 06:02:21 by jwe]
line wrap: on
line diff
--- a/libcruft/ranlib/Basegen.doc Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/Basegen.doc Thu Oct 15 06:02:21 1998 +0000 @@ -9,7 +9,7 @@ - RANLIB + RANDLIB Library of Fortran Routines for Random Number Generation
--- a/libcruft/ranlib/Makefile.in Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/Makefile.in Thu Oct 15 06:02:21 1998 +0000 @@ -12,8 +12,8 @@ top_srcdir = @top_srcdir@ VPATH = @srcdir@ -SPECIAL = README ranlib.chs ranlib.fdoc tstbot.for tstgmn.for \ - tstmid.for +SPECIAL = HOWTOGET README randlib.chs randlib.fdoc \ + tstbot.for tstgmn.for tstmid.for EXTERNAL_DISTFILES = $(DISTFILES) Basegen.doc
--- a/libcruft/ranlib/README Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/README Thu Oct 15 06:02:21 1998 +0000 @@ -9,13 +9,12 @@ - RANLIB + RANDLIB Library of Fortran Routines for Random Number Generation - - + Version 1.3 -- August, 1997 @@ -33,7 +32,8 @@ Barry W. Brown James Lovato - + Kathy Russell + John Venier @@ -52,7 +52,17 @@ This work was supported by grant CA-16672 from the National Cancer Institute. - SUMMARY OF RANLIB + + THANKS TO OUR SUPPORTERS + +This work was supported in part by grant CA-16672 from the National +Cancer Institute. We are grateful to Larry and Pat McNeil of Corpus +Cristi for their generous support. Some equipment used in this effort +was provided by IBM as part of a cooperative study agreement; we thank +them. + + + SUMMARY OF RANDLIB The bottom level routines provide 32 virtual random number generators. Each generator can provide 1,048,576 blocks of numbers, and each block @@ -75,19 +85,23 @@ (10) Random permutations of an integer array (11) Real uniform random deviates between specified limits (12) Binomial random deviates - (13) Poisson random deviates - (14) Integer uniform deviates between specified limits - (15) Seeds for the random number generator calculated from a + (13) Negative Binomial random deviates + (14) Multinomial random deviates + (15) Poisson random deviates + (16) Integer uniform deviates between specified limits + (17) Seeds for the random number generator calculated from a character string INSTALLATION -Directory src contains the Fortran source for most of the routines. -Directory linpack contains two linpack routines needed by the -multivariate generator. If linpack is present on your machine, you -won't need these routines. The Fortran code from these directories -should be compiled and placed in a library. Directory test contains -three test programs for this code. +Directory src contains the Fortran source. The Fortran code from this +directory should be compiled and placed in a library. Directory test +contains three test programs for this code. + + + + + DOCUMENTATION @@ -95,20 +109,19 @@ documentation is in the form of character (ASCII) files. An explanation of the concepts involved in the base generator and details of its implementation are contained in Basegen.doc. A summary of all -of the available routines is contained in ranlib.chs (chs is an +of the available routines is contained in randlib.chs (chs is an abbreviation of 'cheat sheet'). The 'chs' file will probably be the -reference to ranlib that is primarily used. The file, ranlib.fdoc, +reference to randlib that is primarily used. The file, randlib.fdoc, contains all comments heading each routine. There is somewhat more information in 'fdoc' than 'chs', but the additional information consists primarily of references to the literature. - SOURCES The following routines, which were written by others and lightly -modified for consistency in packaging, are included in RANLIB. +modified for consistency in packaging, are included in RANDLIB. Bottom Level Routines @@ -140,6 +153,11 @@ Ahrens, J.H. and Dieter, U. Computer Methods for Sampling from Gamma, Beta, Poisson and Binomial Distributions. Computing, 12 (1974), 223-246. Adaptation of algorithm GS. + + + + + Normal @@ -156,6 +174,7 @@ Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random Variate Generation. Communications of the ACM, 31, 2 (February, 1988) 216. + Poisson This code was obtained from netlib. @@ -163,9 +182,6 @@ Ahrens, J.H. and Dieter, U. Computer Generation of Poisson Deviates From Modified Normal Distributions. ACM Trans. Math. Software, 8, 2 (June 1982),163-179 - - - Beta @@ -183,6 +199,8 @@ Dongarra, J. J., Moler, C. B., Bunch, J. R. and Stewart, G. W. Linpack User's Guide. SIAM Press, Philadelphia. (1979) + + LEGALITIES @@ -203,7 +221,7 @@ Krogh, F. Algorithms Policy. ACM Tran. Math. Softw. 13(1987), 183-186. -We place the Ranlib code that we have written in the public domain. +We place the Randlib code that we have written in the public domain. NO WARRANTY @@ -223,3 +241,106 @@ PARTIES) THE PROGRAM. (Above NO WARRANTY modified from the GNU NO WARRANTY statement.) + + + + WHAT'S NEW IN VERSION 1.1? + + +Random number generation for the Negative Binomial and Multinomial +distributions has been included. + +Two errors in the code which generates random numbers from the Gamma +distribution were fixed. + + + WHAT'S NEW IN VERSION 1.2? + +We changed the name of the package from 'ranlib' to 'randlib'. This +was done so that we can determine who archives it. 'ranlib' is the +name of a Unix utility which produces many spurious hits on a web +search engine. + + +The linpack routines are now housed in the /src directory. + +In several routines, some variables were given an explicit SAVE +attribute and some dummy initial values were changed to prevent +potential errors. +'genbet.f' 'ignbin.f' 'ignpoi.f' 'phrtsd.f' 'sexpo.f' 'sgamma.f' +'snorm.f' + +In several routines, argument checking was implemented; the code now +breaks if inappropriate values are passed to it. +'genbet.f' A and B must be >= 1.0E-37 instead of 0.0 +'genexp.f' AV must be >= 0.0 +'gengam.f' A and R both must be > 0.0 +'gennor.f' SD must be >= 0.0 +'ignbin.f' N must be >= 0, and 0.0 <= PP <= 1.0. +'ignnbn.f' N must be > 0, 0.0 < P < 1.0 (previously allowed N = 0) +'ignpoi.f' MU must be >= 0.0 + +For the Non-Central Chi-Squared and Non-Central F distributions, the +case DF = 1.0 (DFN = 1.0 for the F) is now allowed. +'gennch.f' 'gennf.f' + +Wherever possible, the user-accessible code now calls the base +generators directly. This means improved performance and fewer +dependencies, but the routines should work exactly as before from the +user's point of view. +'genchi.f' 'genf.f' 'gennch.f' 'gennf.f' 'ignnbn.f' + +Many minor modifications have been made which should make the code +more robust, without changing how the code is used. +'genbet.f' 'gengam.f' 'ignpoi.f' 'ignuin.f' 'sgamma.f' 'tstmid.f' + +Finally, five distributions have been added to the mid-level tester, +which test the Exponential, Gamma, Multinomial, Negative Binomial, and +Normal distributions. +'tstmid.f' + + + + + WHAT'S NOT NEW IN VERSION 1.2 ? + +No calling sequences have changed. + + WHAT'S NEW IN VERSION 1.3? + +The calling sequence of SETGMN has been changed! We added an argument +(INTEGER LDCOVM) representing the leading actual dimension of COVM, to +allow the user to use this routine in the case that COVM is contained +in a larger array. This change also makes the routine more compatible +with LINPACK routines. See the following files for details: +'setgmn.f' in the /src directory, and 'randlib.fdoc' and 'randlib.chs' +in the /doc directory. + +Briefly, the declaration of SETGMN has been changed +from: + SUBROUTINE setgmn(meanv,covm,p,parm) +to: + SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm) + +The program 'tstgmn.f' (in the /test directory) was changed to reflect +the change in the calling sequence of SETGMN. + +'randlib.fdoc' and 'randlib.chs' in the /doc directory were changed to +relect the change in the calling sequence of SETGMN. + +Minor changes were made in two routines ('sgamma.f' and 'sexpo.f') to +fix unusual bugs. + +The protection from overflow in deviate generation in two routines +('genf.f' and 'gennf.f') was changed to prevent a constant from +underflowing at compile time. + + WHAT'S NOT NEW IN VERSION 1.3 ? + +No calling sequences (other than SETGMN) have changed. + + MANY THANKS + +The authors would like to thank the many users who have reported bugs +and suggested improvements; Randlib would not be the same today +without them. We heartily encourage others to join them.
--- a/libcruft/ranlib/advnst.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/advnst.f Thu Oct 15 06:02:21 1998 +0000 @@ -60,8 +60,7 @@ IF (qrgnin()) GO TO 10 WRITE (*,*) ' ADVNST called before random number generator ', + ' initialized -- abort!' - CALL XSTOPX - + (' ADVNST called before random number generator initialized') + STOP ' ADVNST called before random number generator initialized' 10 CALL getcgn(g) C
--- a/libcruft/ranlib/genbet.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/genbet.f Thu Oct 15 06:02:21 1998 +0000 @@ -18,16 +18,18 @@ C C A --> First parameter of the beta distribution C REAL A +C JJV (A > 1.0E-37) C C B --> Second parameter of the beta distribution C REAL B +C JJV (B > 1.0E-37) C C C Method C C C R. C. H. Cheng -C Generating Beta Variatew with Nonintegral Shape Parameters +C Generating Beta Variates with Nonintegral Shape Parameters C Communications of the ACM, 21:317-322 (1978) C (Algorithms BB and BC) C @@ -35,10 +37,15 @@ C .. Parameters .. C Close to the largest number that can be exponentiated REAL expmax - PARAMETER (expmax=89.0) +C JJV changed this - 89 was too high, and LOG(1.0E38) = 87.49823 + PARAMETER (expmax=87.49823) C Close to the largest representable single precision number REAL infnty PARAMETER (infnty=1.0E38) +C JJV added the parameter minlog +C Close to the smallest number of which a LOG can be taken. + REAL minlog + PARAMETER (minlog=1.0E-37) C .. C .. Scalar Arguments .. REAL aa,bb @@ -56,18 +63,21 @@ INTRINSIC exp,log,max,min,sqrt C .. C .. Save statement .. - SAVE olda,oldb,alpha,beta,gamma,k1,k2 +C JJV added a,b + SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b C .. C .. Data statements .. - DATA olda,oldb/-1,-1/ +C JJV changed these to ridiculous values + DATA olda,oldb/-1.0E37,-1.0E37/ C .. C .. Executable Statements .. qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb) IF (qsame) GO TO 20 - IF (.NOT. (aa.LE.0.0.OR.bb.LE.0.0)) GO TO 10 - WRITE (*,*) ' AA or BB <= 0 in GENBET - Abort!' +C JJV added small minimum for small log problem in calc of W + IF (.NOT. (aa.LT.minlog.OR.bb.LT.minlog)) GO TO 10 + WRITE (*,*) ' AA or BB < ',minlog,' in GENBET - Abort!' WRITE (*,*) ' AA: ',aa,' BB ',bb - CALL XSTOPX (' AA or BB <= 0 in GENBET - Abort!') + STOP ' AA or BB too small in GENBET - Abort!' 10 olda = aa oldb = bb @@ -92,11 +102,16 @@ C u2 = ranf() v = beta*log(u1/ (1.0-u1)) - IF (.NOT. (v.GT.expmax)) GO TO 50 - w = infnty +C JJV altered this + IF (v.GT.expmax) GO TO 55 +C JJV added checker to see if a*exp(v) will overflow +C JJV 50 _was_ w = a*exp(v); also note here a > 1.0 + 50 w = exp(v) + IF (w.GT.infnty/a) GO TO 55 + w = a*w GO TO 60 + 55 w = infnty - 50 w = a*exp(v) 60 z = u1**2*u2 r = gamma*v - 1.3862944 s = a + r - w @@ -112,6 +127,13 @@ C C Step 4 C +C JJV added checker to see if log(alpha/(b+w)) will +C JJV overflow. If so, we count the log as -INF, and +C JJV consequently evaluate conditional as true, i.e. +C JJV the algorithm rejects the trial and starts over +C JJV May not need this here since ALPHA > 2.0 + IF (alpha/(b+w).LT.minlog) GO TO 40 + IF ((r+alpha*log(alpha/ (b+w))).LT.t) GO TO 40 C C Step 5 @@ -157,12 +179,28 @@ 130 z = u1**2*u2 IF (.NOT. (z.LE.0.25)) GO TO 160 v = beta*log(u1/ (1.0-u1)) - IF (.NOT. (v.GT.expmax)) GO TO 140 - w = infnty - GO TO 150 + +C JJV instead of checking v > expmax at top, I will check +C JJV if a < 1, then check the appropriate values - 140 w = a*exp(v) - 150 GO TO 200 + IF (a.GT.1.0) GO TO 135 +C JJV A < 1 so it can help out if EXP(V) would overflow + IF (v.GT.expmax) GO TO 132 + w = a*exp(v) + GO TO 200 + 132 w = v + log(a) + IF (w.GT.expmax) GO TO 140 + w = exp(w) + GO TO 200 + +C JJV in this case A > 1 + 135 IF (v.GT.expmax) GO TO 140 + w = exp(v) + IF (w.GT.infnty/a) GO TO 140 + w = a*w + GO TO 200 + 140 w = infnty + GO TO 200 160 IF (z.GE.k2) GO TO 120 C @@ -172,12 +210,31 @@ C Step 5 C 170 v = beta*log(u1/ (1.0-u1)) - IF (.NOT. (v.GT.expmax)) GO TO 180 - w = infnty + +C JJV same kind of checking as above + IF (a.GT.1.0) GO TO 175 +C JJV A < 1 so it can help out if EXP(V) would overflow + IF (v.GT.expmax) GO TO 172 + w = a*exp(v) + GO TO 190 + 172 w = v + log(a) + IF (w.GT.expmax) GO TO 180 + w = exp(w) GO TO 190 - 180 w = a*exp(v) - 190 IF ((alpha* (log(alpha/ (b+w))+v)-1.3862944).LT.log(z)) GO TO 120 +C JJV in this case A > 1 + 175 IF (v.GT.expmax) GO TO 180 + w = exp(v) + IF (w.GT.infnty/a) GO TO 180 + w = a*w + GO TO 190 + + 180 w = infnty + +C JJV here we also check to see if log overlows; if so, we treat it +C JJV as -INF, which means condition is true, i.e. restart + 190 IF (alpha/(b+w).LT.minlog) GO TO 120 + IF ((alpha* (log(alpha/ (b+w))+v)-1.3862944).LT.log(z)) GO TO 120 C C Step 6 C
--- a/libcruft/ranlib/genchi.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/genchi.f Thu Oct 15 06:02:21 1998 +0000 @@ -30,16 +30,20 @@ REAL df C .. C .. External Functions .. - REAL gengam - EXTERNAL gengam +C REAL gengam +C EXTERNAL gengam + REAL sgamma + EXTERNAL sgamma C .. C .. Executable Statements .. IF (.NOT. (df.LE.0.0)) GO TO 10 WRITE (*,*) 'DF <= 0 in GENCHI - ABORT' WRITE (*,*) 'Value of DF: ',df - CALL XSTOPX ('DF <= 0 in GENCHI - ABORT') + STOP 'DF <= 0 in GENCHI - ABORT' - 10 genchi = 2.0*gengam(1.0,df/2.0) +C JJV changed this to call sgamma directly +C 10 genchi = 2.0*gengam(1.0,df/2.0) + 10 genchi = 2.0*sgamma(df/2.0) RETURN END
--- a/libcruft/ranlib/genexp.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/genexp.f Thu Oct 15 06:02:21 1998 +0000 @@ -1,4 +1,5 @@ REAL FUNCTION genexp(av) + C********************************************************************** C C REAL FUNCTION GENEXP( AV ) @@ -19,6 +20,7 @@ C AV --> The mean of the exponential distribution from which C a random deviate is to be generated. C REAL AV +C JJV (AV >= 0) C C GENEXP <-- The random deviate. C REAL GENEXP @@ -46,7 +48,13 @@ EXTERNAL sexpo C .. C .. Executable Statements .. - genexp = sexpo()*av +C JJV added check to ensure AV >= 0.0 + IF (av.GE.0.0) GO TO 10 + WRITE (*,*) 'AV < 0.0 in GENEXP - ABORT' + WRITE (*,*) 'Value of AV: ',av + STOP 'AV < 0.0 in GENEXP - ABORT' + + 10 genexp = sexpo()*av RETURN END
--- a/libcruft/ranlib/genf.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/genf.f Thu Oct 15 06:02:21 1998 +0000 @@ -36,24 +36,34 @@ C .. Local Scalars .. REAL xden,xnum C .. +C JJV changed this code to call sgamma directly C .. External Functions .. - REAL genchi - EXTERNAL genchi +C REAL genchi +C EXTERNAL genchi + REAL sgamma + EXTERNAL sgamma C .. C .. Executable Statements .. IF (.NOT. (dfn.LE.0.0.OR.dfd.LE.0.0)) GO TO 10 WRITE (*,*) 'Degrees of freedom nonpositive in GENF - abort!' WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd - CALL XSTOPX ('Degrees of freedom nonpositive in GENF - abort!') + STOP 'Degrees of freedom nonpositive in GENF - abort!' + + 10 xnum = 2.0*sgamma(dfn/2.0)/dfn - 10 xnum = genchi(dfn)/dfn C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) - xden = genchi(dfd)/dfd - IF (.NOT. (xden.LE. (1.2E-38*xnum))) GO TO 20 + xden = 2.0*sgamma(dfd/2.0)/dfd +C JJV changed constant so that it will not underflow at compile time +C JJV while not slowing generator by using double precision or logs. +C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 20 + IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 20 WRITE (*,*) ' GENF - generated numbers would cause overflow' WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden - WRITE (*,*) ' GENF returning 1.0E38' - genf = 1.0E38 +C JJV next 2 lines changed to maintain truncation of large deviates. +C WRITE (*,*) ' GENF returning 1.0E38' +C genf = 1.0E38 + WRITE (*,*) ' GENF returning 1.0E37' + genf = 1.0E37 GO TO 30 20 genf = xnum/xden
--- a/libcruft/ranlib/gengam.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/gengam.f Thu Oct 15 06:02:21 1998 +0000 @@ -16,11 +16,12 @@ C Arguments C C +C JJV added the argument ranges supported C A --> Location parameter of Gamma distribution -C REAL A +C REAL A ( A > 0 ) C C R --> Shape parameter of Gamma distribution -C REAL R +C REAL R ( R > 0 ) C C C Method @@ -37,7 +38,8 @@ C Comm. ACM, 25,1 (Jan. 1982), 47 - 54. C Algorithm GD C -C (Case 0.0 <= R <= 1.0) +C JJV altered the following to reflect sgamma argument ranges +C (Case 0.0 < R < 1.0) C Ahrens, J.H. and Dieter, U. C Computer Methods for Sampling from Gamma, C Beta, Poisson and Binomial Distributions. @@ -53,8 +55,17 @@ EXTERNAL sgamma C .. C .. Executable Statements .. - gengam = sgamma(r) - gengam = gengam/a + +C JJV added argument value checker + IF ( a.GT.0.0 .AND. r.GT.0.0 ) GO TO 10 + WRITE (*,*) 'In GENGAM - Either (1) Location param A <= 0.0 or' + WRITE (*,*) '(2) Shape param R <= 0.0 - ABORT!' + WRITE (*,*) 'A value: ',a,'R value: ',r + STOP 'Location or shape param out of range in GENGAM - ABORT!' +C JJV end addition + + 10 gengam = sgamma(r)/a +C gengam = gengam/a RETURN END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ranlib/genmul.f Thu Oct 15 06:02:21 1998 +0000 @@ -0,0 +1,92 @@ + SUBROUTINE genmul(n,p,ncat,ix) +C********************************************************************** +C +C SUBROUTINE GENMUL( N, P, NCAT, IX ) +C GENerate an observation from the MULtinomial distribution +C +C +C Arguments +C +C +C N --> Number of events that will be classified into one of +C the categories 1..NCAT +C INTEGER N +C +C P --> Vector of probabilities. P(i) is the probability that +C an event will be classified into category i. Thus, P(i) +C must be [0,1]. Only the first NCAT-1 P(i) must be defined +C since P(NCAT) is 1.0 minus the sum of the first +C NCAT-1 P(i). +C REAL P(NCAT-1) +C +C NCAT --> Number of categories. Length of P and IX. +C INTEGER NCAT +C +C IX <-- Observation from multinomial distribution. All IX(i) +C will be nonnegative and their sum will be N. +C INTEGER IX(NCAT) +C +C +C Method +C +C +C Algorithm from page 559 of +C +C Devroye, Luc +C +C Non-Uniform Random Variate Generation. Springer-Verlag, +C New York, 1986. +C +C********************************************************************** +C .. Scalar Arguments .. + INTEGER n,ncat +C .. +C .. Array Arguments .. + REAL p(*) + INTEGER ix(*) +C .. +C .. Local Scalars .. + REAL prob,ptot,sum + INTEGER i,icat,ntot +C .. +C .. External Functions .. + INTEGER ignbin + EXTERNAL ignbin +C .. +C .. Intrinsic Functions .. + INTRINSIC abs +C .. +C .. Executable Statements .. + +C Check Arguments + IF (n.LT.0) STOP 'N < 0 in GENMUL' + IF (ncat.LE.1) STOP 'NCAT <= 1 in GENMUL' + ptot = 0.0 + DO 10,i = 1,ncat - 1 + IF (p(i).LT.0.0) STOP 'Some P(i) < 0 in GENMUL' + IF (p(i).GT.1.0) STOP 'Some P(i) > 1 in GENMUL' + ptot = ptot + p(i) + 10 CONTINUE + IF (ptot.GT.0.99999) STOP 'Sum of P(i) > 1 in GENMUL' + +C Initialize variables + ntot = n + sum = 1.0 + DO 20,i = 1,ncat + ix(i) = 0 + 20 CONTINUE + +C Generate the observation + DO 30,icat = 1,ncat - 1 + prob = p(icat)/sum + ix(icat) = ignbin(ntot,prob) + ntot = ntot - ix(icat) + IF (ntot.LE.0) RETURN + sum = sum - p(icat) + 30 CONTINUE + ix(ncat) = ntot + +C Finished + RETURN + + END
--- a/libcruft/ranlib/gennch.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/gennch.f Thu Oct 15 06:02:21 1998 +0000 @@ -18,7 +18,7 @@ C C C DF --> Degrees of freedom of the chisquare -C (Must be > 1.0) +C (Must be >= 1.0) C REAL DF C C XNONC --> Noncentrality parameter of the chisquare @@ -31,26 +31,39 @@ C C Uses fact that noncentral chisquare is the sum of a chisquare C deviate with DF-1 degrees of freedom plus the square of a normal -C deviate with mean XNONC and standard deviation 1. +C deviate with mean sqrt(XNONC) and standard deviation 1. C C********************************************************************** C .. Scalar Arguments .. REAL df,xnonc C .. C .. External Functions .. - REAL genchi,gennor - EXTERNAL genchi,gennor +C JJV changed these to call SGAMMA and SNORM directly +C REAL genchi,gennor +C EXTERNAL genchi,gennor + REAL sgamma,snorm + EXTERNAL sgamma,snorm C .. C .. Intrinsic Functions .. INTRINSIC sqrt C .. +C JJV changed abort to df < 1, and added case: df = 1 C .. Executable Statements .. - IF (.NOT. (df.LE.1.0.OR.xnonc.LT.0.0)) GO TO 10 - WRITE (*,*) 'DF <= 1 or XNONC < 0 in GENNCH - ABORT' + IF (.NOT. (df.LT.1.0.OR.xnonc.LT.0.0)) GO TO 10 + WRITE (*,*) 'DF < 1 or XNONC < 0 in GENNCH - ABORT' WRITE (*,*) 'Value of DF: ',df,' Value of XNONC',xnonc - CALL XSTOPX ('DF <= 1 or XNONC < 0 in GENNCH - ABORT') + STOP 'DF < 1 or XNONC < 0 in GENNCH - ABORT' + +C JJV changed this to call SGAMMA and SNORM directly +C gennch = genchi(df-1.0) + gennor(sqrt(xnonc),1.0)**2 - 10 gennch = genchi(df-1.0) + gennor(sqrt(xnonc),1.0)**2 - RETURN + 10 IF (df.GE.1.000001) GO TO 20 +C JJV case DF = 1.0 + gennch = (snorm() + sqrt(xnonc))**2 + GO TO 30 +C JJV case DF > 1.0 + 20 gennch = 2.0*sgamma((df-1.0)/2.0) + (snorm() + sqrt(xnonc))**2 + 30 RETURN + END
--- a/libcruft/ranlib/gennf.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/gennf.f Thu Oct 15 06:02:21 1998 +0000 @@ -45,31 +45,52 @@ LOGICAL qcond C .. C .. External Functions .. - REAL genchi,gennch - EXTERNAL genchi,gennch +C JJV changed the code to call SGAMMA and SNORM directly +C REAL genchi,gennch +C EXTERNAL genchi,gennch + REAL sgamma,snorm + EXTERNAL sgamma,snorm C .. C .. Executable Statements .. - qcond = dfn .LE. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0 +C JJV changed the argument checker to allow DFN = 1.0 +C JJV in the same way as GENNCH was changed. + qcond = dfn .LT. 1.0 .OR. dfd .LE. 0.0 .OR. xnonc .LT. 0.0 IF (.NOT. (qcond)) GO TO 10 - WRITE (*,*) 'In GENNF - Either (1) Numerator DF <= 1.0 or' - WRITE (*,*) '(2) Denominator DF < 0.0 or ' + WRITE (*,*) 'In GENNF - Either (1) Numerator DF < 1.0 or' + WRITE (*,*) '(2) Denominator DF <= 0.0 or ' WRITE (*,*) '(3) Noncentrality parameter < 0.0' WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd,'XNONC value: ', + xnonc - CALL XSTOPX - + ('Degrees of freedom or noncent param our of range in GENNF') + STOP 'Degrees of freedom or noncent param out of range in GENNF' + +C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) +C JJV changed this to call SGAMMA and SNORM directly +C xnum = gennch(dfn,xnonc)/dfn + 10 IF (dfn.GE.1.000001) GO TO 20 +C JJV case dfn = 1.0 - here I am treating dfn as exactly 1.0 + xnum = (snorm() + sqrt(xnonc))**2 + GO TO 30 - 10 xnum = gennch(dfn,xnonc)/dfn -C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) - xden = genchi(dfd)/dfd - IF (.NOT. (xden.LE. (1.2E-38*xnum))) GO TO 20 +C JJV case dfn > 1.0 + 20 xnum = (2.0*sgamma((dfn-1.0)/2.0) + (snorm()+sqrt(xnonc))**2)/dfn + +C xden = genchi(dfd)/dfd + 30 xden = 2.0*sgamma(dfd/2.0)/dfd + +C JJV changed constant so that it will not underflow at compile time +C JJV while not slowing generator by using double precision or logs. +C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 40 + IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 40 WRITE (*,*) ' GENNF - generated numbers would cause overflow' WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden - WRITE (*,*) ' GENNF returning 1.0E38' - gennf = 1.0E38 - GO TO 30 +C JJV next 2 lines changed to maintain truncation of large deviates. +C WRITE (*,*) ' GENNF returning 1.0E38' +C gennf = 1.0E38 + WRITE (*,*) ' GENNF returning 1.0E37' + gennf = 1.0E37 + GO TO 50 - 20 gennf = xnum/xden - 30 RETURN + 40 gennf = xnum/xden + 50 RETURN END
--- a/libcruft/ranlib/gennor.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/gennor.f Thu Oct 15 06:02:21 1998 +0000 @@ -21,6 +21,7 @@ C C SD --> Standard deviation of the normal distribution. C REAL SD +C JJV (SD >= 0) C C GENNOR <-- Generated normal deviate. C REAL GENNOR @@ -48,7 +49,13 @@ EXTERNAL snorm C .. C .. Executable Statements .. - gennor = sd*snorm() + av +C JJV added check to ensure SD >= 0.0 + IF (sd.GE.0.0) GO TO 10 + WRITE (*,*) 'SD < 0.0 in GENNOR - ABORT' + WRITE (*,*) 'Value of SD: ',sd + STOP 'SD < 0.0 in GENNOR - ABORT' + + 10 gennor = sd*snorm() + av RETURN END
--- a/libcruft/ranlib/genunf.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/genunf.f Thu Oct 15 06:02:21 1998 +0000 @@ -33,7 +33,7 @@ IF (.NOT. (low.GT.high)) GO TO 10 WRITE (*,*) 'LOW > HIGH in GENUNF: LOW ',low,' HIGH: ',high WRITE (*,*) 'Abort' - CALL XSTOPX ('LOW > High in GENUNF - Abort') + STOP 'LOW > High in GENUNF - Abort' 10 genunf = low + (high-low)*ranf()
--- a/libcruft/ranlib/getcgn.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/getcgn.f Thu Oct 15 06:02:21 1998 +0000 @@ -47,7 +47,7 @@ IF (.NOT. (g.LT.0.OR.g.GT.numg)) GO TO 10 WRITE (*,*) ' Generator number out of range in SETCGN:', + ' Legal range is 1 to ',numg,' -- ABORT!' - CALL XSTOPX (' Generator number out of range in SETCGN') + STOP ' Generator number out of range in SETCGN' 10 curntg = g RETURN
--- a/libcruft/ranlib/getsd.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/getsd.f Thu Oct 15 06:02:21 1998 +0000 @@ -1,7 +1,7 @@ SUBROUTINE getsd(iseed1,iseed2) C********************************************************************** C -C SUBROUTINE GETSD(ISEED1,ISEED2) +C SUBROUTINE GETSD(G,ISEED1,ISEED2) C GET SeeD C C Returns the value of two integer seeds of the current generator @@ -62,8 +62,7 @@ IF (qrgnin()) GO TO 10 WRITE (*,*) ' GETSD called before random number generator ', + ' initialized -- abort!' - CALL XSTOPX - + (' GETSD called before random number generator initialized') + STOP ' GETSD called before random number generator initialized' 10 CALL getcgn(g) iseed1 = cg1(g)
--- a/libcruft/ranlib/ignbin.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/ignbin.f Thu Oct 15 06:02:21 1998 +0000 @@ -1,7 +1,7 @@ INTEGER FUNCTION ignbin(n,pp) C********************************************************************** C -C INTEGER FUNCTION IGNBIN( N, P ) +C INTEGER FUNCTION IGNBIN( N, PP ) C C GENerate BINomial random deviate C @@ -20,11 +20,13 @@ C N --> The number of trials in the binomial distribution C from which a random deviate is to be generated. C INTEGER N +C JJV (N >= 0) C -C P --> The probability of an event in each trial of the +C PP --> The probability of an event in each trial of the C binomial distribution from which a random deviate C is to be generated. -C REAL P +C REAL PP +C JJV (0.0 <= pp <= 1.0) C C IGNBIN <-- A random deviate yielding the number of events C from N independent trials, each of which has @@ -162,9 +164,16 @@ C .. C .. Intrinsic Functions .. INTRINSIC abs,alog,amin1,iabs,int,sqrt +C JJV .. +C JJV .. Save statement .. + SAVE p,q,m,fm,xnp,xnpq,p1,xm,xl,xr,c,xll,xlr,p2,p3,p4,qn,r,g, + + psave,nsave +C JJV I am including the variables in data statements C .. C .. Data statements .. - DATA psave,nsave/-1.,-1/ +C JJV made these ridiculous starting values - the hope is that +C JJV no one will call this the first time with them as args + DATA psave,nsave/-1.0E37,-214748365/ C .. C .. Executable Statements .. IF (pp.NE.psave) GO TO 10 @@ -173,10 +182,18 @@ C C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE C - 10 psave = pp + +C JJV added the argument checker - involved only renaming 10 +C JJV and 20 to the checkers and adding checkers +C JJV Only remaining problem - if called initially with the +C JJV initial values of psave and nsave, it will hang + 10 IF (pp.LT.0.0) STOP 'PP < 0.0 in IGNBIN - ABORT!' + IF (pp.GT.1.0) STOP 'PP > 1.0 in IGNBIN - ABORT!' + psave = pp p = amin1(psave,1.-psave) q = 1. - p - 20 xnp = n*p + 20 IF (n.LT.0) STOP 'N < 0 in IGNBIN - ABORT!' + xnp = n*p nsave = n IF (xnp.LT.30.) GO TO 140 ffm = xnp + p
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ranlib/ignnbn.f Thu Oct 15 06:02:21 1998 +0000 @@ -0,0 +1,78 @@ + INTEGER FUNCTION ignnbn(n,p) +C********************************************************************** +C +C INTEGER FUNCTION IGNNBN( N, P ) +C +C GENerate Negative BiNomial random deviate +C +C +C Function +C +C +C Generates a single random deviate from a negative binomial +C distribution. +C +C +C Arguments +C +C +C N --> Required number of events. +C INTEGER N +C JJV (N > 0) +C +C P --> The probability of an event during a Bernoulli trial. +C REAL P +C JJV (0.0 < P < 1.0) +C +C +C +C Method +C +C +C Algorithm from page 480 of +C +C Devroye, Luc +C +C Non-Uniform Random Variate Generation. Springer-Verlag, +C New York, 1986. +C +C********************************************************************** +C .. +C .. Scalar Arguments .. + REAL p + INTEGER n +C .. +C .. Local Scalars .. + REAL y,a,r +C .. +C .. External Functions .. +C JJV changed to call SGAMMA directly +C REAL gengam + REAL sgamma + INTEGER ignpoi +C EXTERNAL gengam,ignpoi + EXTERNAL sgamma,ignpoi +C .. +C .. Intrinsic Functions .. + INTRINSIC real +C .. +C .. Executable Statements .. +C Check Arguments +C JJV changed argumnet checker to abort if N <= 0 + IF (n.LE.0) STOP 'N <= 0 in IGNNBN' + IF (p.LE.0.0) STOP 'P <= 0.0 in IGNNBN' + IF (p.GE.1.0) STOP 'P >= 1.0 in IGNNBN' + +C Generate Y, a random gamma (n,(1-p)/p) variable +C JJV Note: the above parametrization is consistent with Devroye, +C JJV but gamma (p/(1-p),n) is the equivalent in our code + 10 r = real(n) + a = p/ (1.0-p) +C y = gengam(a,r) + y = sgamma(r)/a + +C Generate a random Poisson(y) variable + ignnbn = ignpoi(y) + RETURN + + END
--- a/libcruft/ranlib/ignpoi.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/ignpoi.f Thu Oct 15 06:02:21 1998 +0000 @@ -1,7 +1,7 @@ INTEGER FUNCTION ignpoi(mu) C********************************************************************** C -C INTEGER FUNCTION IGNPOI( AV ) +C INTEGER FUNCTION IGNPOI( MU ) C C GENerate POIsson random deviate C @@ -10,18 +10,19 @@ C C C Generates a single random deviate from a Poisson -C distribution with mean AV. +C distribution with mean MU. C C C Arguments C C -C AV --> The mean of the Poisson distribution from which +C MU --> The mean of the Poisson distribution from which C a random deviate is to be generated. -C REAL AV +C REAL MU +C JJV (MU >= 0.0) C -C GENEXP <-- The random deviate. -C REAL GENEXP +C IGNPOI <-- The random deviate. +C INTEGER IGNPOI (non-negative) C C C Method @@ -68,7 +69,7 @@ C C C -C MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B. +C MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B C TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT C COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL C @@ -82,7 +83,8 @@ C .. Local Scalars .. REAL a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e, + fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx - INTEGER j,k,kflag,l,m +C JJV I added a variable 'll' here - it is the 'l' for CASE A + INTEGER j,k,kflag,l,ll,m C .. C .. Local Arrays .. REAL fact(10),pp(35) @@ -94,27 +96,40 @@ C .. Intrinsic Functions .. INTRINSIC abs,alog,exp,float,ifix,max0,min0,sign,sqrt C .. +C JJV added this for case: mu unchanged +C .. Save statement .. + SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0, + + a0, a1, a2, a3, a4, a5, a6, a7, fact, muprev, muold +C .. +C JJV end addition - I am including vars in Data statements C .. Data statements .. - DATA muprev,muold/0.,0./ +C JJV changed initial values of MUPREV and MUOLD to -1.0E37 +C JJV if no one calls IGNPOI with MU = -1.0E37 the first time, +C JJV the code shouldn't break + DATA muprev,muold/-1.0E37,-1.0E37/ DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118, + -.1661269,.1421878,-.1384794,.1250060/ DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./ C .. C .. Executable Statements .. + IF (mu.EQ.muprev) GO TO 10 IF (mu.LT.10.0) GO TO 120 C -C C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED) +C C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED) C +C JJV This is the case where I changed 'l' to 'll' +C JJV Here 'll' is set once and used in a comparison once + muprev = mu s = sqrt(mu) d = 6.0*mu*mu C C THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL -C PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484) +C PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484) C IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . C - l = ifix(mu-1.1484) + ll = ifix(mu-1.1484) C C STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE C @@ -124,7 +139,7 @@ C C STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH C - IF (ignpoi.GE.l) RETURN + IF (ignpoi.GE.ll) RETURN C C STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U C @@ -214,9 +229,16 @@ C C C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) C - 120 muprev = 0.0 +C JJV changed MUPREV assignment from 0.0 to initial value + 120 muprev = -1.0E37 IF (mu.EQ.muold) GO TO 130 - muold = mu +C JJV added argument checker here + IF (mu.GE.0.0) GO TO 125 + WRITE (*,*) 'MU < 0 in IGNPOI - ABORT' + WRITE (*,*) 'Value of MU: ',mu + STOP 'MU < 0 in IGNPOI - ABORT' +C JJV added line label here + 125 muold = mu m = max0(1,ifix(mu)) l = 0 p = exp(-mu)
--- a/libcruft/ranlib/ignuin.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/ignuin.f Thu Oct 15 06:02:21 1998 +0000 @@ -68,8 +68,6 @@ ignuin = low RETURN - GO TO 70 - C Number to be generated should be in range 0..RANGE C Set MAXNOW so that the number of integers in 0..MAXNOW is an C integral multiple of the number in 0..RANGE @@ -77,14 +75,10 @@ 30 ranp1 = range + 1 maxnow = (maxnum/ranp1)*ranp1 40 ign = ignlgi() - 1 - IF (.NOT. (ign.LE.maxnow)) GO TO 50 + IF (.NOT. (ign.LE.maxnow)) GO TO 40 ignuin = low + mod(ign,ranp1) RETURN - 50 GO TO 40 - - 60 CONTINUE - 70 CONTINUE 80 IF (.NOT. (err.EQ.1)) GO TO 90 WRITE (*,*) err1 GO TO 100 @@ -94,10 +88,7 @@ 100 WRITE (*,*) ' LOW: ',low,' HIGH: ',high WRITE (*,*) ' Abort on Fatal ERROR' IF (.NOT. (err.EQ.1)) GO TO 110 - CALL XSTOPX ('LOW > HIGH in IGNUIN') - IGNUIN = 0 - - GO TO 120 + STOP 'LOW > HIGH in IGNUIN' 110 STOP ' ( HIGH - LOW ) > 2,147,483,561 in IGNUIN'
--- a/libcruft/ranlib/initgn.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/initgn.f Thu Oct 15 06:02:21 1998 +0000 @@ -66,8 +66,7 @@ IF (qrgnin()) GO TO 10 WRITE (*,*) ' INITGN called before random number generator ', + ' initialized -- abort!' - CALL XSTOPX - + (' INITGN called before random number generator initialized') + STOP ' INITGN called before random number generator initialized' 10 CALL getcgn(g) IF ((-1).NE. (isdtyp)) GO TO 20
--- a/libcruft/ranlib/mltmod.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/mltmod.f Thu Oct 15 06:02:21 1998 +0000 @@ -39,7 +39,7 @@ WRITE (*,*) ' A, M, S out of order in MLTMOD - ABORT!' WRITE (*,*) ' A = ',a,' S = ',s,' M = ',m WRITE (*,*) ' MLTMOD requires: 0 < A < M; 0 < S < M' - CALL XSTOPX (' A, M, S out of order in MLTMOD - ABORT!') + STOP ' A, M, S out of order in MLTMOD - ABORT!' 10 IF (.NOT. (a.LT.h)) GO TO 20 a0 = a
--- a/libcruft/ranlib/phrtsd.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/phrtsd.f Thu Oct 15 06:02:21 1998 +0000 @@ -59,6 +59,11 @@ C .. Intrinsic Functions .. INTRINSIC index,mod C .. +C JJV added Save statement for variable in Data statement +C .. Save statements .. + SAVE shift +C JJV end addition +C .. C .. Data statements .. DATA shift/1,64,4096,262144,16777216/ C ..
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ranlib/randlib.chs Thu Oct 15 06:02:21 1998 +0000 @@ -0,0 +1,368 @@ + SUMMARY OF ROUTINES IN RANDLIB + +0. Base Level Routines to Set and Obtain Values of Seeds + +(These should be the only base level routines used by those who don't +need multiple generators with blocks of numbers.) + +C********************************************************************** +C +C SUBROUTINE SETALL(ISEED1,ISEED2) +C SET ALL random number generators +C INTEGER ISEED1, ISEED2 +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GETSD(ISEED1,ISEED2) +C GET SeeD +C INTEGER ISEED1, ISEED2 +C +C Returns the value of two integer seeds of the current generator +C in ISEED1, ISEED2 +C +C********************************************************************** + +I. Higher Level Routines + +C********************************************************************** +C +C REAL FUNCTION GENBET( A, B ) +C GeNerate BETa random deviate +C REAL A,B +C +C Returns a single random deviate from the beta distribution with +C parameters A and B. The density of the beta is +C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENCHI( DF ) +C Generate random value of CHIsquare variable +C REAL DF +C +C Generates random deviate from the distribution of a chisquare +C with DF degrees of freedom random variable. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENEXP( AV ) +C GENerate EXPonential random deviate +C REAL AV +C +C Generates a single random deviate from an exponential +C distribution with mean AV. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENF( DFN, DFD ) +C GENerate random deviate from the F distribution +C REAL DFN, DFD +C +C Generates a random deviate from the F (variance ratio) +C distribution with DFN degrees of freedom in the numerator +C and DFD degrees of freedom in the denominator. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENGAM( A, R ) +C GENerates random deviates from GAMma distribution +C REAL A, R +C +C Generates random deviates from the gamma distribution whose +C density is +C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENMN(PARM,X,WORK) +C GENerate Multivariate Normal random deviate +C REAL PARM(*), X(*), WORK(*) +C +C PARM is set by SETGMN which must be called prior to GENMN. The +C generated deviates are placed in X. WORK is a work array of the +C same size as X. +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENMUL( N, P, NCAT, IX ) +C GENerate MULtinomial random deviate +C REAL P(*) +C INTEGER N, NCAT, IX(*) +C +C Generates deviates from a Multinomial distribution with NCAT +C categories. P specifies the probability of an event in each +C category. The generated deviates are placed in IX. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNCH( DF, XNONC ) +C Generate random value of Noncentral CHIsquare variable +C REAL DF, XNONC +C +C Generates random deviate from the distribution of a noncentral +C chisquare with DF degrees of freedom and noncentrality parameter +C XNONC. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNF( DFN, DFD, XNONC ) +C GENerate random deviate from the Noncentral F distribution +C REAL DFN, DFD, XNONC +C +C Generates a random deviate from the noncentral F (variance ratio) +C distribution with DFN degrees of freedom in the numerator, and DFD +C degrees of freedom in the denominator, and noncentrality parameter +C XNONC. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNOR( AV, SD ) +C GENerate random deviate from a NORmal distribution +C REAL AV, SD +C +C Generates a single random deviate from a normal distribution +C with mean, AV, and standard deviation, SD. +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENPRM( IARRAY, LARRAY ) +C GENerate random PeRMutation of iarray +C INTEGER IARRAY(LARRAY), LARRAY +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENUNF( LOW, HIGH ) +C GeNerate Uniform Real between LOW and HIGH +C REAL LOW, HIGH +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNBIN( N, P ) +C GENerate BINomial random deviate +C INTEGER N +C REAL P +C +C Returns a single random deviate from a binomial +C distribution whose number of trials is N and whose +C probability of an event in each trial is P. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNNBN( N, P ) +C GENerate Negative BiNomial random deviate +C INTEGER N +C REAL P +C +C Returns a single random deviate from a negative binomial +C distribution with number of events N and whose +C probability of an event in each trial is P. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNPOI( AV ) +C GENerate POIsson random deviate +C REAL AV +C +C Generates a single random deviate from a Poisson +C distribution with mean AV. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNUIN( LOW, HIGH ) +C GeNerate Uniform INteger +C INTEGER LOW, HIGH +C +C Generates an integer uniformly distributed between LOW and HIGH. +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE PHRTSD( PHRASE, SEED1, SEED2 ) +C PHRase To SeeDs +C CHARACTER*(*) PHRASE +C INTEGER SEED1, SEED2 +C +C Uses a phrase (character string) to generate two seeds for the RGN +C random number generator. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION RANF() +C RANDom number generator as a Function +C +C Returns a random floating point number from a uniform distribution +C over 0 - 1 (endpoints of this interval are not returned) using the +C current generator +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM) +C SET Generate Multivariate Normal random deviate +C INTEGER LDCOVM, P +C REAL MEANV(P), COVM(LDCOVM,P), PARM(P*(P+3)/2 + 1) +C +C P is the length of normal vectors to be generated, MEANV +C is the vector of their means and COVM(1:P,1:P) is their variance +C covariance matrix. LDCOVM is the leading actual dimension of +C COVM, which this routine needs to know although only the +C (1:P,1:P) slice of COVM is used. +C Places information necessary to generate the deviates in PARM. +C +C********************************************************************** + +II. Uniform Generator and Associated Routines + + + A. SETTING THE SEED OF ALL GENERATORS + +C********************************************************************** +C +C SUBROUTINE SETALL(ISEED1,ISEED2) +C SET ALL random number generators +C INTEGER ISEED1, ISEED2 +C +C********************************************************************** + + B. OBTAINING RANDOM NUMBERS + +C********************************************************************** +C +C INTEGER FUNCTION IGNLGI() +C GeNerate LarGe Integer +C +C Returns a random integer following a uniform distribution over +C (1, 2147483562) using the current generator. +C +C********************************************************************** + +C********************************************************************** +C +C REAL FUNCTION RANF() +C RANDom number generator as a Function +C +C Returns a random floating point number from a uniform distribution +C over 0 - 1 (endpoints of this interval are not returned) using the +C current generator +C +C********************************************************************** + + C. SETTING AND OBTAINING THE NUMBER OF THE CURRENT GENERATOR + +C********************************************************************** +C +C SUBROUTINE SETCGN( G ) +C Set GeNerator +C INTEGER G +C +C Sets the current generator to G. All references to a generator +C are to the current generator. +C +C********************************************************************** + +C********************************************************************** +C +C SUBROUTINE GETCGN(G) +C GET Current GeNerator +C INTEGER G +C +C Returns in G the number of the current random number generator +C +C********************************************************************** + + D. OBTAINING OR CHANGING SEEDS IN CURRENT GENERATOR + +C********************************************************************** +C +C SUBROUTINE ADVNST(K) +C ADV-a-N-ce ST-ate +C INTEGER K +C +C Advances the state of the current generator by 2^K values and +C resets the initial seed to that value. +C +C********************************************************************** + +C********************************************************************** +C +C SUBROUTINE GETSD(ISEED1,ISEED2) +C GET SeeD +C INTEGER ISEED1, ISEED2 +C +C Returns the value of two integer seeds of the current generator +C in ISEED1, ISEED2 +C +C********************************************************************** + +C********************************************************************** +C +C SUBROUTINE INITGN(ISDTYP) +C INIT-ialize current G-e-N-erator +C +C INTEGER ISDTYP The state to which the generator is to be set +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********************************************************************** + +C********************************************************************** +C +C SUBROUTINE SETSD(ISEED1,ISEED2) +C SET S-ee-D of current generator +C +C Resets the initial seed of the current generator to ISEED1 and +C ISEED2. The seeds of the other generators remain unchanged. +C +C********************************************************************** + + E. MISCELLANY + +C********************************************************************** +C +C INTEGER FUNCTION MLTMOD(A,S,M) +C Returns (A*S) MOD M +C INTEGER A, S, M +C +C********************************************************************** + +C********************************************************************** +C +C SUBROUTINE SETANT(QVALUE) +C SET ANTithetic +C LOGICAL QVALUE +C +C Sets whether the current generator produces antithetic values. If +C X is the value normally returned from a uniform [0,1] random +C number generator then 1 - X is the antithetic value. If X is the +C value normally returned from a uniform [0,N] random number +C generator then N - 1 - X is the antithetic value. +C +C All generators are initialized to NOT generate antithetic values. +C +C********************************************************************** + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ranlib/randlib.fdoc Thu Oct 15 06:02:21 1998 +0000 @@ -0,0 +1,961 @@ + + + + + + + + + + + + RANDLIB + + Library of Fortran Routines for Random Number Generation + + + + + + + + + Full Documentation of Each Routine + + + + + + + + + Compiled and Written by: + + Barry W. Brown + James Lovato + + + + + + + + + + + Department of Biomathematics, Box 237 + The University of Texas, M.D. Anderson Cancer Center + 1515 Holcombe Boulevard + Houston, TX 77030 + + + This work was supported by grant CA-16672 from the National Cancer Institute. + +C********************************************************************** +C +C SUBROUTINE ADVNST(K) +C ADV-a-N-ce ST-ate +C +C Advances the state of the current generator by 2^K values and +C resets the initial seed to that value. +C +C This is a transcription from Pascal to Fortran of routine +C Advance_State 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 K -> The generator is advanced by2^K values +C INTEGER K +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENBET( A, B ) +C GeNerate BETa random deviate +C +C +C Function +C +C +C Returns a single random deviate from the beta distribution with +C parameters A and B. The density of the beta is +C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 +C +C +C Arguments +C +C +C A --> First parameter of the beta distribution +C REAL A +C (A >= 1.0E-37) +C +C B --> Second parameter of the beta distribution +C REAL B +C (B >= 1.0E-37) +C +C +C Method +C +C +C R. C. H. Cheng +C Generating Beta Variables with Nonintegral Shape Parameters +C Communications of the ACM, 21:317-322 (1978) +C (Algorithms BB and BC) +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENCHI( DF ) +C Generate random value of CHIsquare variable +C +C +C Function +C +C +C Generates random deviate from the distribution of a chisquare +C with DF degrees of freedom random variable. +C +C +C Arguments +C +C +C DF --> Degrees of freedom of the chisquare +C (Must be positive) +C REAL DF +C +C +C Method +C +C +C Uses relation between chisquare and gamma. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENEXP( AV ) +C +C GENerate EXPonential random deviate +C +C +C Function +C +C +C Generates a single random deviate from an exponential +C distribution with mean AV. +C +C +C Arguments +C +C +C AV --> The mean of the exponential distribution from which +C a random deviate is to be generated. +C REAL AV +C (AV >= 0) +C +C GENEXP <-- The random deviate. +C REAL GENEXP +C +C +C Method +C +C +C Renames SEXPO from TOMS as slightly modified by BWB to use RANF +C instead of SUNIF. +C +C For details see: +C +C Ahrens, J.H. and Dieter, U. +C Computer Methods for Sampling From the +C Exponential and Normal Distributions. +C Comm. ACM, 15,10 (Oct. 1972), 873 - 882. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENF( DFN, DFD ) +C GENerate random deviate from the F distribution +C +C +C Function +C +C +C Generates a random deviate from the F (variance ratio) +C distribution with DFN degrees of freedom in the numerator +C and DFD degrees of freedom in the denominator. +C +C +C Arguments +C +C +C DFN --> Numerator degrees of freedom +C (Must be positive) +C REAL DFN +C DFD --> Denominator degrees of freedom +C (Must be positive) +C REAL DFD +C +C +C Method +C +C +C Directly generates ratio of chisquare variates +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENGAM( A, R ) +C GENerates random deviates from GAMma distribution +C +C +C Function +C +C +C Generates random deviates from the gamma distribution whose +C density is +C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) +C +C +C Arguments +C +C +C A --> Location parameter of Gamma distribution +C REAL A ( A > 0 ) +C +C R --> Shape parameter of Gamma distribution +C REAL R ( R > 0 ) +C +C +C Method +C +C +C Renames SGAMMA from TOMS as slightly modified by BWB to use RANF +C instead of SUNIF. +C +C For details see: +C (Case R >= 1.0) +C Ahrens, J.H. and Dieter, U. +C Generating Gamma Variates by a +C Modified Rejection Technique. +C Comm. ACM, 25,1 (Jan. 1982), 47 - 54. +C Algorithm GD +C +C (Case 0.0 < R < 1.0) +C Ahrens, J.H. and Dieter, U. +C Computer Methods for Sampling from Gamma, +C Beta, Poisson and Binomial Distributions. +C Computing, 12 (1974), 223-246/ +C Adapted algorithm GS. +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENMN(PARM,X,WORK) +C GENerate Multivariate Normal random deviate +C +C +C Arguments +C +C +C PARM --> Parameters needed to generate multivariate normal +C deviates (MEANV and Cholesky decomposition of +C COVM). Set by a previous call to SETGMN. +C +C 1 : 1 - size of deviate, P +C 2 : P + 1 - mean vector +C P+2 : P*(P+3)/2 + 1 - upper half of cholesky +C decomposition of cov matrix +C REAL PARM(*) +C +C X <-- Vector deviate generated. +C REAL X(P) +C +C WORK <--> Scratch array +C REAL WORK(P) +C +C +C Method +C +C +C 1) Generate P independent standard normal deviates - Ei ~ N(0,1) +C +C 2) SETGMN uses Cholesky decomposition find A s.t. trans(A)*A = COV +C +C 3) Generate trans(A)*E + MEANV ~ N(MEANV,COVM) +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENMUL( N, P, NCAT, IX ) +C GENerate an observation from the MULtinomial distribution +C +C +C Arguments +C +C +C N --> Number of events that will be classified into one of +C the categories 1..NCAT +C INTEGER N +C (N >= 0) +C +C P --> Vector of probabilities. P(i) is the probability that +C an event will be classified into category i. Thus, P(i) +C must be [0,1]. Only the first NCAT-1 P(i) must be defined +C since P(NCAT) is 1.0 minus the sum of the first +C NCAT-1 P(i). +C REAL P(NCAT-1) +C +C NCAT --> Number of categories. Length of P and IX. +C INTEGER NCAT +C (NCAT > 1) +C +C IX <-- Observation from multinomial distribution. All IX(i) +C will be nonnegative and their sum will be N. +C INTEGER IX(NCAT) +C +C +C Method +C +C +C Algorithm from page 559 of +C +C Devroye, Luc +C +C Non-Uniform Random Variate Generation. Springer-Verlag, +C New York, 1986. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNCH( DF, XNONC ) +C Generate random value of Noncentral CHIsquare variable +C +C +C Function +C +C +C +C Generates random deviate from the distribution of a noncentral +C chisquare with DF degrees of freedom and noncentrality parameter +C XNONC. +C +C +C Arguments +C +C +C DF --> Degrees of freedom of the chisquare +C (Must be >= 1.0) +C REAL DF +C +C XNONC --> Noncentrality parameter of the chisquare +C (Must be >= 0.0) +C REAL XNONC +C +C +C Method +C +C +C Uses fact that noncentral chisquare is the sum of a chisquare +C deviate with DF-1 degrees of freedom plus the square of a normal +C deviate with mean XNONC and standard deviation 1. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNF( DFN, DFD, XNONC ) +C GENerate random deviate from the Noncentral F distribution +C +C +C Function +C +C +C Generates a random deviate from the noncentral F (variance ratio) +C distribution with DFN degrees of freedom in the numerator, and DFD +C degrees of freedom in the denominator, and noncentrality parameter +C XNONC. +C +C +C Arguments +C +C +C DFN --> Numerator degrees of freedom +C (Must be >= 1.0) +C REAL DFN +C DFD --> Denominator degrees of freedom +C (Must be positive) +C REAL DFD +C +C XNONC --> Noncentrality parameter +C (Must be nonnegative) +C REAL XNONC +C +C +C Method +C +C +C Directly generates ratio of noncentral numerator chisquare variate +C to central denominator chisquare variate. +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENNOR( AV, SD ) +C +C GENerate random deviate from a NORmal distribution +C +C +C Function +C +C +C Generates a single random deviate from a normal distribution +C with mean, AV, and standard deviation, SD. +C +C +C Arguments +C +C +C AV --> Mean of the normal distribution. +C REAL AV +C +C SD --> Standard deviation of the normal distribution. +C REAL SD +C (SD >= 0) +C +C GENNOR <-- Generated normal deviate. +C REAL GENNOR +C +C +C Method +C +C +C Renames SNORM from TOMS as slightly modified by BWB to use RANF +C instead of SUNIF. +C +C For details see: +C Ahrens, J.H. and Dieter, U. +C Extensions of Forsythe's Method for Random +C Sampling from the Normal Distribution. +C Math. Comput., 27,124 (Oct. 1973), 927 - 937. +C +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GENPRM( IARRAY, LARRAY ) +C GENerate random PeRMutation of iarray +C +C +C Arguments +C +C +C IARRAY <--> On output IARRAY is a random permutation of its +C value on input +C INTEGER IARRAY( LARRAY ) +C +C LARRAY <--> Length of IARRAY +C INTEGER LARRAY +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION GENUNF( LOW, HIGH ) +C +C GeNerate Uniform Real between LOW and HIGH +C +C +C Function +C +C +C Generates a real uniformly distributed between LOW and HIGH. +C +C +C Arguments +C +C +C LOW --> Low bound (exclusive) on real value to be generated +C REAL LOW +C +C HIGH --> High bound (exclusive) on real value to be generated +C REAL HIGH +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GETCGN(G) +C Get GeNerator +C +C Returns in G the number of the current random number generator +C +C +C Arguments +C +C +C G <-- Number of the current random number generator (1..32) +C INTEGER G +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE GETSD(ISEED1,ISEED2) +C GET SeeD +C +C Returns the value of two integer seeds of the current generator +C +C This is a transcription from Pascal to Fortran of routine +C Get_State 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 +C ISEED1 <- First integer seed of generator G +C INTEGER ISEED1 +C +C ISEED2 <- Second integer seed of generator G +C INTEGER ISEED1 +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNBIN( N, P ) +C +C GENerate BINomial random deviate +C +C +C Function +C +C +C Generates a single random deviate from a binomial +C distribution whose number of trials is N and whose +C probability of an event in each trial is P. +C +C +C Arguments +C +C +C N --> The number of trials in the binomial distribution +C from which a random deviate is to be generated. +C INTEGER N +C (N >= 0) +C +C P --> The probability of an event in each trial of the +C binomial distribution from which a random deviate +C is to be generated. +C REAL P +C (0.0 <= P <= 1.0) +C +C IGNBIN <-- A random deviate yielding the number of events +C from N independent trials, each of which has +C a probability of event P. +C INTEGER IGNBIN +C +C +C Note +C +C +C Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set +C by a call similar to the following +C DUM = RANSET( ISEED1, ISEED2 ) +C +C +C Method +C +C +C This is algorithm BTPE from: +C +C Kachitvichyanukul, V. and Schmeiser, B. W. +C +C Binomial Random Variate Generation. +C Communications of the ACM, 31, 2 +C (February, 1988) 216. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNNBN( N, P ) +C +C GENerate Negative BiNomial random deviate +C +C +C Function +C +C +C Generates a single random deviate from a negative binomial +C distribution. +C +C +C Arguments +C +C +C N --> Required number of events. +C INTEGER N +C (N > 0) +C +C P --> The probability of an event during a Bernoulli trial. +C REAL P +C (0.0 < P < 1.0) +C +C +C +C Method +C +C +C Algorithm from page 480 of +C +C Devroye, Luc +C +C Non-Uniform Random Variate Generation. Springer-Verlag, +C New York, 1986. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNLGI() +C GeNerate LarGe Integer +C +C Returns a random integer following a uniform distribution over +C (1, 2147483562) using the current generator. +C +C This is a transcription from Pascal to Fortran of routine +C Random 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********************************************************************** +C +C INTEGER FUNCTION IGNPOI( MU ) +C +C GENerate POIsson random deviate +C +C +C Function +C +C +C Generates a single random deviate from a Poisson +C distribution with mean MU. +C +C +C Arguments +C +C +C MU --> The mean of the Poisson distribution from which +C a random deviate is to be generated. +C REAL MU +C (MU >= 0.0) +C +C IGNPOI <-- The random deviate. +C REAL IGNPOI (non-negative) +C +C +C Method +C +C +C Renames KPOIS from TOMS as slightly modified by BWB to use RANF +C instead of SUNIF. +C +C For details see: +C +C Ahrens, J.H. and Dieter, U. +C Computer Generation of Poisson Deviates +C From Modified Normal Distributions. +C ACM Trans. Math. Software, 8, 2 +C (June 1982),163-179 +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION IGNUIN( LOW, HIGH ) +C +C GeNerate Uniform INteger +C +C +C Function +C +C +C Generates an integer uniformly distributed between LOW and HIGH. +C +C +C Arguments +C +C +C LOW --> Low bound (inclusive) on integer value to be generated +C INTEGER LOW +C +C HIGH --> High bound (inclusive) on integer value to be generated +C INTEGER HIGH +C +C +C Note +C +C +C If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and +C stops the program. +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE INITGN(ISDTYP) +C INIT-ialize current G-e-N-erator +C +C Reinitializes the state of the current generator +C ISDTYP = -1 => sets the state to its initial seed +C ISDTYP = 0 => sets the state to its last (previous) seed +C ISDTYP = 1 => sets the state to a new seed 2^w values +C from its last seed +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 INTEGER ISDTYP +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE INRGCM() +C INitialize Random number Generator CoMmon +C +C +C Function +C +C +C Initializes common area for random number generator. This saves +C the nuisance of a BLOCK DATA routine and the difficulty of +C assuring that the routine is loaded with the other routines. +C +C********************************************************************** +C********************************************************************** +C +C INTEGER FUNCTION MLTMOD(A,S,M) +C +C Returns (A*S) MOD M +C +C This is a transcription from Pascal to Fortran of routine +C MULtMod_Decompos 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 A, S, M --> +C INTEGER A,S,M +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE PHRTSD( PHRASE, SEED1, SEED2 ) +C PHRase To SeeDs +C +C +C Function +C +C +C Uses a phrase (character string) to generate two seeds for the RGN +C random number generator. +C +C +C Arguments +C +C +C PHRASE --> Phrase to be used for random number generation +C CHARACTER*(*) PHRASE +C +C SEED1 <-- First seed for RGN generator +C INTEGER SEED1 +C +C SEED2 <-- Second seed for RGN generator +C INTEGER SEED2 +C +C +C Note +C +C +C Trailing blanks are eliminated before the seeds are generated. +C +C Generated seed values will fall in the range 1..2^30 +C (1..1,073,741,824) +C +C********************************************************************** +C********************************************************************** +C +C REAL FUNCTION RANF() +C RANDom number generator as a Function +C +C Returns a random floating point number from a uniform distribution +C over 0 - 1 (endpoints of this interval are not returned) using the +C current generator +C +C This is a transcription from Pascal to Fortran of routine +C Uniform_01 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********************************************************************** +C +C SUBROUTINE SETALL(ISEED1,ISEED2) +C SET ALL random number generators +C +C Sets the initial seed of generator 1 to ISEED1 and ISEED2. The +C initial seeds of the other generators are set accordingly, and +C all generators states are set to these seeds. +C +C This is a transcription from Pascal to Fortran of routine +C Set_Initial_Seed 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 ISEED1 -> First of two integer seeds +C INTEGER ISEED1 +C +C ISEED2 -> Second of two integer seeds +C INTEGER ISEED1 +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE SETANT(QVALUE) +C SET ANTithetic +C +C Sets whether the current generator produces antithetic values. If +C X is the value normally returned from a uniform [0,1] random +C number generator then 1 - X is the antithetic value. If X is the +C value normally returned from a uniform [0,N] random number +C generator then N - 1 - X is the antithetic value. +C +C All generators are initialized to NOT generate antithetic values. +C +C This is a transcription from Pascal to Fortran of routine +C Set_Antithetic 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 QVALUE -> .TRUE. if generator G is to generating antithetic +C values, otherwise .FALSE. +C LOGICAL QVALUE +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE SETCGN( G ) +C Set GeNerator +C +C Sets the current generator to G. All references to a generato +C are to the current generator. +C +C +C Arguments +C +C +C G --> Number of the current random number generator (1..32) +C INTEGER G +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM) +C SET Generate Multivariate Normal random deviate +C +C +C Function +C +C +C Places P, MEANV, and the Cholesky factoriztion of COVM +C in PARM for GENMN. +C +C +C Arguments +C +C +C MEANV --> Mean vector of multivariate normal distribution. +C REAL MEANV(P) +C +C COVM <--> (Input) Covariance matrix of the multivariate +C normal distribution. This routine uses only the +C (1:P,1:P) slice of COVM, but needs to know LDCOVM. +C +C (Output) Destroyed on output +C REAL COVM(LDCOVM,P) +C +C LDCOVM --> Leading actual dimension of COVM. +C INTEGER LDCOVM +C +C P --> Dimension of the normal, or length of MEANV. +C INTEGER P +C +C PARM <-- Array of parameters needed to generate multivariate +C normal deviates (P, MEANV and Cholesky decomposition +C of COVM). +C 1 : 1 - P +C 2 : P + 1 - MEANV +C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM +C REAL PARM(P*(P+3)/2 + 1) +C +C********************************************************************** +C********************************************************************** +C +C SUBROUTINE SETSD(ISEED1,ISEED2) +C SET S-ee-D of current generator +C +C Resets the initial seed and state of generator g to ISEED1 and +C ISEED2. The seeds and states of the other generators remain +C unchanged. +C +C This is a transcription from Pascal to Fortran of routine +C Set_Seed 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 ISEED1 -> First integer seed +C INTEGER ISEED1 +C +C ISEED2 -> Second integer seed +C INTEGER ISEED1 +C +C**********************************************************************
--- a/libcruft/ranlib/ranlib.chs Fri Oct 09 03:24:56 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,336 +0,0 @@ - SUMMARY OF ROUTINES IN RANLIB - -0. Base Level Routines to Set and Obtain Values of Seeds - -(These should be the only base level routines used by those who don't -need multiple generators with blocks of numbers.) - -C********************************************************************** -C -C SUBROUTINE SETALL(ISEED1,ISEED2) -C SET ALL random number generators -C INTEGER ISEED1, ISEED2 -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GETSD(ISEED1,ISEED2) -C GET SeeD -C INTEGER ISEED1, ISEED2 -C -C Returns the value of two integer seeds of the current generator -C in ISEED1, ISEED2 -C -C********************************************************************** - -I. Higher Level Routines - -C********************************************************************** -C -C REAL FUNCTION GENBET( A, B ) -C GeNerate BETa random deviate -C REAL A,B -C -C Returns a single random deviate from the beta distribution with -C parameters A and B. The density of the beta is -C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENCHI( DF ) -C Generate random value of CHIsquare variable -C REAL DF -C -C Generates random deviate from the distribution of a chisquare -C with DF degrees of freedom random variable. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENEXP( AV ) -C GENerate EXPonential random deviate -C REAL AV -C -C Generates a single random deviate from an exponential -C distribution with mean AV. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENF( DFN, DFD ) -C GENerate random deviate from the F distribution -C REAL DFN, DFD -C -C Generates a random deviate from the F (variance ratio) -C distribution with DFN degrees of freedom in the numerator -C and DFD degrees of freedom in the denominator. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENGAM( A, R ) -C GENerates random deviates from GAMma distribution -C REAL A, R -C -C Generates random deviates from the gamma distribution whose -C density is -C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GENMN(PARM,X,WORK) -C GENerate Multivariate Normal random deviate -C REAL PARM(*), X(*), WORK(*) -C -C PARM is set by SETGMN which must be called prior to GENMN. The -C generated deviates are placed in X. WORK is a work array of the -C same size as X. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNCH( DF, XNONC ) -C Generate random value of Noncentral CHIsquare variable -C REAL DF, XNONC -C -C Generates random deviate from the distribution of a noncentral -C chisquare with DF degrees of freedom and noncentrality parameter -C XNONC. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNF( DFN, DFD, XNONC ) -C GENerate random deviate from the Noncentral F distribution -C REAL DFN, DFD, XNONC -C -C Generates a random deviate from the noncentral F (variance ratio) -C distribution with DFN degrees of freedom in the numerator, and DFD -C degrees of freedom in the denominator, and noncentrality parameter -C XNONC. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNOR( AV, SD ) -C GENerate random deviate from a NORmal distribution -C REAL AV, SD -C -C Generates a single random deviate from a normal distribution -C with mean, AV, and standard deviation, SD. -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GENPRM( IARRAY, LARRAY ) -C GENerate random PeRMutation of iarray -C INTEGER IARRAY(LARRAY), LARRAY -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENUNF( LOW, HIGH ) -C GeNerate Uniform Real between LOW and HIGH -C REAL LOW, HIGH -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNBIN( N, P ) -C GENerate BINomial random deviate -C INTEGER N -C REAL P -C -C Returns a single random deviate from a binomial -C distribution whose number of trials is N and whose -C probability of an event in each trial is P. -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNPOI( AV ) -C GENerate POIsson random deviate -C REAL AV -C -C Generates a single random deviate from a Poisson -C distribution with mean AV. -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNUIN( LOW, HIGH ) -C GeNerate Uniform INteger -C INTEGER LOW, HIGH -C -C Generates an integer uniformly distributed between LOW and HIGH. -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE PHRTSD( PHRASE, SEED1, SEED2 ) -C PHRase To SeeDs -C CHARACTER*(*) PHRASE -C INTEGER SEED1, SEED2 -C -C Uses a phrase (character string) to generate two seeds for the RGN -C random number generator. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION RANF() -C RANDom number generator as a Function -C -C Returns a random floating point number from a uniform distribution -C over 0 - 1 (endpoints of this interval are not returned) using the -C current generator -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE SETGMN( MEANV, COVM, P, PARM) -C SET Generate Multivariate Normal random deviate -C INTEGER P -C REAL MEANV(P), COVM(P,P), PARM(P*(P+3)/2 + 1) -C -C P is the length of normal vectors to be generated, MEANV -C is the vector of their means and COVM is their variance -C covariance matrix. Places information necessary to generate -C the deviates in PARM. -C -C********************************************************************** - -II. Uniform Generator and Associated Routines - - - A. SETTING THE SEED OF ALL GENERATORS - -C********************************************************************** -C -C SUBROUTINE SETALL(ISEED1,ISEED2) -C SET ALL random number generators -C INTEGER ISEED1, ISEED2 -C -C********************************************************************** - - B. OBTAINING RANDOM NUMBERS - -C********************************************************************** -C -C INTEGER FUNCTION IGNLGI() -C GeNerate LarGe Integer -C -C Returns a random integer following a uniform distribution over -C (1, 2147483562) using the current generator. -C -C********************************************************************** - -C********************************************************************** -C -C REAL FUNCTION RANF() -C RANDom number generator as a Function -C -C Returns a random floating point number from a uniform distribution -C over 0 - 1 (endpoints of this interval are not returned) using the -C current generator -C -C********************************************************************** - - C. SETTING AND OBTAINING THE NUMBER OF THE CURRENT GENERATOR - -C********************************************************************** -C -C SUBROUTINE SETCGN( G ) -C Set GeNerator -C INTEGER G -C -C Sets the current generator to G. All references to a generator -C are to the current generator. -C -C********************************************************************** - -C********************************************************************** -C -C SUBROUTINE GETCGN(G) -C GET Current GeNerator -C INTEGER G -C -C Returns in G the number of the current random number generator -C -C********************************************************************** - - D. OBTAINING OR CHANGING SEEDS IN CURRENT GENERATOR - -C********************************************************************** -C -C SUBROUTINE ADVNST(K) -C ADV-a-N-ce ST-ate -C INTEGER K -C -C Advances the state of the current generator by 2^K values and -C resets the initial seed to that value. -C -C********************************************************************** - -C********************************************************************** -C -C SUBROUTINE GETSD(ISEED1,ISEED2) -C GET SeeD -C INTEGER ISEED1, ISEED2 -C -C Returns the value of two integer seeds of the current generator -C in ISEED1, ISEED2 -C -C********************************************************************** - -C********************************************************************** -C -C SUBROUTINE INITGN(ISDTYP) -C INIT-ialize current G-e-N-erator -C -C INTEGER ISDTYP The state to which the generator is to be set -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********************************************************************** - -C********************************************************************** -C -C SUBROUTINE SETSD(ISEED1,ISEED2) -C SET S-ee-D of current generator -C -C Resets the initial seed of the current generator to ISEED1 and -C ISEED2. The seeds of the other generators remain unchanged. -C -C********************************************************************** - - E. MISCELLANY - -C********************************************************************** -C -C INTEGER FUNCTION MLTMOD(A,S,M) -C Returns (A*S) MOD M -C INTEGER A, S, M -C -C********************************************************************** - -C********************************************************************** -C -C SUBROUTINE SETANT(QVALUE) -C SET ANTithetic -C LOGICAL QVALUE -C -C Sets whether the current generator produces antithetic values. If -C X is the value normally returned from a uniform [0,1] random -C number generator then 1 - X is the antithetic value. If X is the -C value normally returned from a uniform [0,N] random number -C generator then N - 1 - X is the antithetic value. -C -C All generators are initialized to NOT generate antithetic values. -C -C**********************************************************************
--- a/libcruft/ranlib/ranlib.fdoc Fri Oct 09 03:24:56 1998 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,870 +0,0 @@ - - - - - - - - - - - - RANLIB - - Library of Fortran Routines for Random Number Generation - - - - - - - - - Full Documentation of Each Routine - - - - - - - - - Compiled and Written by: - - Barry W. Brown - James Lovato - - - - - - - - - - - Department of Biomathematics, Box 237 - The University of Texas, M.D. Anderson Cancer Center - 1515 Holcombe Boulevard - Houston, TX 77030 - - - This work was supported by grant CA-16672 from the National Cancer Institute. - -C********************************************************************** -C -C SUBROUTINE ADVNST(K) -C ADV-a-N-ce ST-ate -C -C Advances the state of the current generator by 2^K values and -C resets the initial seed to that value. -C -C This is a transcription from Pascal to Fortran of routine -C Advance_State 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 K -> The generator is advanced by2^K values -C INTEGER K -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENBET( A, B ) -C GeNerate BETa random deviate -C -C -C Function -C -C -C Returns a single random deviate from the beta distribution with -C parameters A and B. The density of the beta is -C x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1 -C -C -C Arguments -C -C -C A --> First parameter of the beta distribution -C REAL A -C -C B --> Second parameter of the beta distribution -C REAL B -C -C -C Method -C -C -C R. C. H. Cheng -C Generating Beta Variables with Nonintegral Shape Parameters -C Communications of the ACM, 21:317-322 (1978) -C (Algorithms BB and BC) -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENCHI( DF ) -C Generate random value of CHIsquare variable -C -C -C Function -C -C -C Generates random deviate from the distribution of a chisquare -C with DF degrees of freedom random variable. -C -C -C Arguments -C -C -C DF --> Degrees of freedom of the chisquare -C (Must be positive) -C REAL DF -C -C -C Method -C -C -C Uses relation between chisquare and gamma. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENEXP( AV ) -C -C GENerate EXPonential random deviate -C -C -C Function -C -C -C Generates a single random deviate from an exponential -C distribution with mean AV. -C -C -C Arguments -C -C -C AV --> The mean of the exponential distribution from which -C a random deviate is to be generated. -C REAL AV -C -C GENEXP <-- The random deviate. -C REAL GENEXP -C -C -C Method -C -C -C Renames SEXPO from TOMS as slightly modified by BWB to use RANF -C instead of SUNIF. -C -C For details see: -C -C Ahrens, J.H. and Dieter, U. -C Computer Methods for Sampling From the -C Exponential and Normal Distributions. -C Comm. ACM, 15,10 (Oct. 1972), 873 - 882. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENF( DFN, DFD ) -C GENerate random deviate from the F distribution -C -C -C Function -C -C -C Generates a random deviate from the F (variance ratio) -C distribution with DFN degrees of freedom in the numerator -C and DFD degrees of freedom in the denominator. -C -C -C Arguments -C -C -C DFN --> Numerator degrees of freedom -C (Must be positive) -C REAL DFN -C DFD --> Denominator degrees of freedom -C (Must be positive) -C REAL DFD -C -C -C Method -C -C -C Directly generates ratio of chisquare variates -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENGAM( A, R ) -C GENerates random deviates from GAMma distribution -C -C -C Function -C -C -C Generates random deviates from the gamma distribution whose -C density is -C (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X) -C -C -C Arguments -C -C -C A --> Location parameter of Gamma distribution -C REAL A -C -C R --> Shape parameter of Gamma distribution -C REAL R -C -C -C Method -C -C -C Renames SGAMMA from TOMS as slightly modified by BWB to use RANF -C instead of SUNIF. -C -C For details see: -C (Case R >= 1.0) -C Ahrens, J.H. and Dieter, U. -C Generating Gamma Variates by a -C Modified Rejection Technique. -C Comm. ACM, 25,1 (Jan. 1982), 47 - 54. -C Algorithm GD -C -C (Case 0.0 <= R <= 1.0) -C Ahrens, J.H. and Dieter, U. -C Computer Methods for Sampling from Gamma, -C Beta, Poisson and Binomial Distributions. -C Computing, 12 (1974), 223-246/ -C Adapted algorithm GS. -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GENMN(PARM,X,WORK) -C GENerate Multivariate Normal random deviate -C -C -C Arguments -C -C -C PARM --> Parameters needed to generate multivariate normal -C deviates (MEANV and Cholesky decomposition of -C COVM). Set by a previous call to SETGMN. -C -C 1 : 1 - size of deviate, P -C 2 : P + 1 - mean vector -C P+2 : P*(P+3)/2 + 1 - upper half of cholesky -C decomposition of cov matrix -C REAL PARM(*) -C -C X <-- Vector deviate generated. -C REAL X(P) -C -C WORK <--> Scratch array -C REAL WORK(P) -C -C -C Method -C -C -C 1) Generate P independent standard normal deviates - Ei ~ N(0,1) -C -C 2) SETGMN uses Cholesky decomposition find A s.t. trans(A)*A = COV -C -C 3) Generate trans(A)*E + MEANV ~ N(MEANV,COVM) -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNCH( DF, XNONC ) -C Generate random value of Noncentral CHIsquare variable -C -C -C Function -C -C - -C Generates random deviate from the distribution of a noncentral -C chisquare with DF degrees of freedom and noncentrality parameter -C XNONC. -C -C -C Arguments -C -C -C DF --> Degrees of freedom of the chisquare -C (Must be > 1.0) -C REAL DF -C -C XNONC --> Noncentrality parameter of the chisquare -C (Must be >= 0.0) -C REAL XNONC -C -C -C Method -C -C -C Uses fact that noncentral chisquare is the sum of a chisquare -C deviate with DF-1 degrees of freedom plus the square of a normal -C deviate with mean XNONC and standard deviation 1. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNF( DFN, DFD, XNONC ) -C GENerate random deviate from the Noncentral F distribution -C -C -C Function -C -C -C Generates a random deviate from the noncentral F (variance ratio) -C distribution with DFN degrees of freedom in the numerator, and DFD -C degrees of freedom in the denominator, and noncentrality parameter -C XNONC. -C -C -C Arguments -C -C -C DFN --> Numerator degrees of freedom -C (Must be >= 1.0) -C REAL DFN -C DFD --> Denominator degrees of freedom -C (Must be positive) -C REAL DFD -C -C XNONC --> Noncentrality parameter -C (Must be nonnegative) -C REAL XNONC -C -C -C Method -C -C -C Directly generates ratio of noncentral numerator chisquare variate -C to central denominator chisquare variate. -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENNOR( AV, SD ) -C -C GENerate random deviate from a NORmal distribution -C -C -C Function -C -C -C Generates a single random deviate from a normal distribution -C with mean, AV, and standard deviation, SD. -C -C -C Arguments -C -C -C AV --> Mean of the normal distribution. -C REAL AV -C -C SD --> Standard deviation of the normal distribution. -C REAL SD -C -C GENNOR <-- Generated normal deviate. -C REAL GENNOR -C -C -C Method -C -C -C Renames SNORM from TOMS as slightly modified by BWB to use RANF -C instead of SUNIF. -C -C For details see: -C Ahrens, J.H. and Dieter, U. -C Extensions of Forsythe's Method for Random -C Sampling from the Normal Distribution. -C Math. Comput., 27,124 (Oct. 1973), 927 - 937. -C -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GENPRM( IARRAY, LARRAY ) -C GENerate random PeRMutation of iarray -C -C -C Arguments -C -C -C IARRAY <--> On output IARRAY is a random permutation of its -C value on input -C INTEGER IARRAY( LARRAY ) -C -C LARRAY <--> Length of IARRAY -C INTEGER LARRAY -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION GENUNF( LOW, HIGH ) -C -C GeNerate Uniform Real between LOW and HIGH -C -C -C Function -C -C -C Generates a real uniformly distributed between LOW and HIGH. -C -C -C Arguments -C -C -C LOW --> Low bound (exclusive) on real value to be generated -C REAL LOW -C -C HIGH --> High bound (exclusive) on real value to be generated -C REAL HIGH -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GETCGN(G) -C Get GeNerator -C -C Returns in G the number of the current random number generator -C -C -C Arguments -C -C -C G <-- Number of the current random number generator (1..32) -C INTEGER G -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE GETSD(ISEED1,ISEED2) -C GET SeeD -C -C Returns the value of two integer seeds of the current generator -C -C This is a transcription from Pascal to Fortran of routine -C Get_State 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 -C ISEED1 <- First integer seed of generator G -C INTEGER ISEED1 -C -C ISEED2 <- Second integer seed of generator G -C INTEGER ISEED1 -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNBIN( N, P ) -C -C GENerate BINomial random deviate -C -C -C Function -C -C -C Generates a single random deviate from a binomial -C distribution whose number of trials is N and whose -C probability of an event in each trial is P. -C -C -C Arguments -C -C -C N --> The number of trials in the binomial distribution -C from which a random deviate is to be generated. -C INTEGER N -C -C P --> The probability of an event in each trial of the -C binomial distribution from which a random deviate -C is to be generated. -C REAL P -C -C IGNBIN <-- A random deviate yielding the number of events -C from N independent trials, each of which has -C a probability of event P. -C INTEGER IGNBIN -C -C -C Note -C -C -C Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set -C by a call similar to the following -C DUM = RANSET( ISEED1, ISEED2 ) -C -C -C Method -C -C -C This is algorithm BTPE from: -C -C Kachitvichyanukul, V. and Schmeiser, B. W. -C -C Binomial Random Variate Generation. -C Communications of the ACM, 31, 2 -C (February, 1988) 216. -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNLGI() -C GeNerate LarGe Integer -C -C Returns a random integer following a uniform distribution over -C (1, 2147483562) using the current generator. -C -C This is a transcription from Pascal to Fortran of routine -C Random 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********************************************************************** -C -C INTEGER FUNCTION IGNPOI( AV ) -C -C GENerate POIsson random deviate -C -C -C Function -C -C -C Generates a single random deviate from a Poisson -C distribution with mean AV. -C -C -C Arguments -C -C -C AV --> The mean of the Poisson distribution from which -C a random deviate is to be generated. -C REAL AV -C -C GENEXP <-- The random deviate. -C REAL GENEXP -C -C -C Method -C -C -C Renames KPOIS from TOMS as slightly modified by BWB to use RANF -C instead of SUNIF. -C -C For details see: -C -C Ahrens, J.H. and Dieter, U. -C Computer Generation of Poisson Deviates -C From Modified Normal Distributions. -C ACM Trans. Math. Software, 8, 2 -C (June 1982),163-179 -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION IGNUIN( LOW, HIGH ) -C -C GeNerate Uniform INteger -C -C -C Function -C -C -C Generates an integer uniformly distributed between LOW and HIGH. -C -C -C Arguments -C -C -C LOW --> Low bound (inclusive) on integer value to be generated -C INTEGER LOW -C -C HIGH --> High bound (inclusive) on integer value to be generated -C INTEGER HIGH -C -C -C Note -C -C -C If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and -C stops the program. -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE INITGN(ISDTYP) -C INIT-ialize current G-e-N-erator -C -C Reinitializes the state of the current generator -C ISDTYP = -1 => sets the state to its initial seed -C ISDTYP = 0 => sets the state to its last (previous) seed -C ISDTYP = 1 => sets the state to a new seed 2^w values -C from its last seed -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 INTEGER ISDTYP -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE INRGCM() -C INitialize Random number Generator CoMmon -C -C -C Function -C -C -C Initializes common area for random number generator. This saves -C the nuisance of a BLOCK DATA routine and the difficulty of -C assuring that the routine is loaded with the other routines. -C -C********************************************************************** -C********************************************************************** -C -C INTEGER FUNCTION MLTMOD(A,S,M) -C -C Returns (A*S) MOD M -C -C This is a transcription from Pascal to Fortran of routine -C MULtMod_Decompos 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 A, S, M --> -C INTEGER A,S,M -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE PHRTSD( PHRASE, SEED1, SEED2 ) -C PHRase To SeeDs -C -C -C Function -C -C -C Uses a phrase (character string) to generate two seeds for the RGN -C random number generator. -C -C -C Arguments -C -C -C PHRASE --> Phrase to be used for random number generation -C CHARACTER*(*) PHRASE -C -C SEED1 <-- First seed for RGN generator -C INTEGER SEED1 -C -C SEED2 <-- Second seed for RGN generator -C INTEGER SEED2 -C -C -C Note -C -C -C Trailing blanks are eliminated before the seeds are generated. -C -C Generated seed values will fall in the range 1..2^30 -C (1..1,073,741,824) -C -C********************************************************************** -C********************************************************************** -C -C REAL FUNCTION RANF() -C RANDom number generator as a Function -C -C Returns a random floating point number from a uniform distribution -C over 0 - 1 (endpoints of this interval are not returned) using the -C current generator -C -C This is a transcription from Pascal to Fortran of routine -C Uniform_01 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********************************************************************** -C -C SUBROUTINE SETALL(ISEED1,ISEED2) -C SET ALL random number generators -C -C Sets the initial seed of generator 1 to ISEED1 and ISEED2. The -C initial seeds of the other generators are set accordingly, and -C all generators states are set to these seeds. -C -C This is a transcription from Pascal to Fortran of routine -C Set_Initial_Seed 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 ISEED1 -> First of two integer seeds -C INTEGER ISEED1 -C -C ISEED2 -> Second of two integer seeds -C INTEGER ISEED1 -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE SETANT(QVALUE) -C SET ANTithetic -C -C Sets whether the current generator produces antithetic values. If -C X is the value normally returned from a uniform [0,1] random -C number generator then 1 - X is the antithetic value. If X is the -C value normally returned from a uniform [0,N] random number -C generator then N - 1 - X is the antithetic value. -C -C All generators are initialized to NOT generate antithetic values. -C -C This is a transcription from Pascal to Fortran of routine -C Set_Antithetic 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 QVALUE -> .TRUE. if generator G is to generating antithetic -C values, otherwise .FALSE. -C LOGICAL QVALUE -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE SETCGN( G ) -C Set GeNerator -C -C Sets the current generator to G. All references to a generato -C are to the current generator. -C -C -C Arguments -C -C -C G --> Number of the current random number generator (1..32) -C INTEGER G -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE SETGMN( MEANV, COVM, P, PARM) -C SET Generate Multivariate Normal random deviate -C -C -C Function -C -C -C Places P, MEANV, and the Cholesky factoriztion of COVM -C in GENMN. -C -C -C Arguments -C -C -C MEANV --> Mean vector of multivariate normal distribution. -C REAL MEANV(P) -C -C COVM <--> (Input) Covariance matrix of the multivariate -C normal distribution -C (Output) Destroyed on output -C REAL COVM(P,P) -C -C P --> Dimension of the normal, or length of MEANV. -C INTEGER P -C -C PARM <-- Array of parameters needed to generate multivariate norma -C deviates (P, MEANV and Cholesky decomposition of -C COVM). -C 1 : 1 - P -C 2 : P + 1 - MEANV -C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM -C REAL PARM(P*(P+3)/2 + 1) -C -C********************************************************************** -C********************************************************************** -C -C SUBROUTINE SETSD(ISEED1,ISEED2) -C SET S-ee-D of current generator -C -C Resets the initial seed and state of generator g to ISEED1 and -C ISEED2. The seeds and states of the other generators remain -C unchanged. -C -C This is a transcription from Pascal to Fortran of routine -C Set_Seed 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 ISEED1 -> First integer seed -C INTEGER ISEED1 -C -C ISEED2 -> Second integer seed -C INTEGER ISEED1 -C -C**********************************************************************
--- a/libcruft/ranlib/setant.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/setant.f Thu Oct 15 06:02:21 1998 +0000 @@ -65,8 +65,7 @@ IF (qrgnin()) GO TO 10 WRITE (*,*) ' SETANT called before random number generator ', + ' initialized -- abort!' - CALL XSTOPX - + (' SETANT called before random number generator initialized') + STOP ' SETANT called before random number generator initialized' 10 CALL getcgn(g) qanti(g) = qvalue
--- a/libcruft/ranlib/setgmn.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/setgmn.f Thu Oct 15 06:02:21 1998 +0000 @@ -1,7 +1,13 @@ - SUBROUTINE setgmn(meanv,covm,p,parm) + SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm) +C SUBROUTINE setgmn(meanv,covm,p,parm) +C JJV changed this routine to take leading dimension of COVM +C JJV argument and pass it to SPOFA, making it easier to use +C JJV if the COVM which is used is contained in a larger matrix +C JJV and to make the routine more consistent with LINPACK. +C JJV Changes are in comments, declarations, and the call to SPOFA. C********************************************************************** C -C SUBROUTINE SETGMN( MEANV, COVM, P, PARM) +C SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM) C SET Generate Multivariate Normal random deviate C C @@ -9,7 +15,7 @@ C C C Places P, MEANV, and the Cholesky factoriztion of COVM -C in GENMN. +C in PARM for GENMN. C C C Arguments @@ -19,16 +25,21 @@ C REAL MEANV(P) C C COVM <--> (Input) Covariance matrix of the multivariate -C normal distribution +C normal distribution. This routine uses only the +C (1:P,1:P) slice of COVM, but needs to know LDCOVM. +C C (Output) Destroyed on output -C REAL COVM(P,P) +C REAL COVM(LDCOVM,P) +C +C LDCOVM --> Leading actual dimension of COVM. +C INTEGER LDCOVM C C P --> Dimension of the normal, or length of MEANV. C INTEGER P C -C PARM <-- Array of parameters needed to generate multivariate norma -C deviates (P, MEANV and Cholesky decomposition of -C COVM). +C PARM <-- Array of parameters needed to generate multivariate +C normal deviates (P, MEANV and Cholesky decomposition +C of COVM). C 1 : 1 - P C 2 : P + 1 - MEANV C P+2 : P*(P+3)/2 + 1 - Cholesky decomposition of COVM @@ -36,10 +47,12 @@ C C********************************************************************** C .. Scalar Arguments .. - INTEGER p +C INTEGER p + INTEGER p, ldcovm C .. C .. Array Arguments .. - REAL covm(p,p),meanv(p),parm(p* (p+3)/2+1) +C REAL covm(p,p),meanv(p),parm(p* (p+3)/2+1) + REAL covm(ldcovm,p),meanv(p),parm(p* (p+3)/2+1) C .. C .. Local Scalars .. INTEGER i,icount,info,j @@ -55,7 +68,7 @@ IF (.NOT. (p.LE.0)) GO TO 10 WRITE (*,*) 'P nonpositive in SETGMN' WRITE (*,*) 'Value of P: ',p - CALL XSTOPX ('P nonpositive in SETGMN') + STOP 'P nonpositive in SETGMN' 10 parm(1) = p C @@ -67,10 +80,11 @@ C C Cholesky decomposition to find A s.t. trans(A)*(A) = COVM C - CALL spofa(covm,p,p,info) +C CALL spofa(covm,p,p,info) + CALL spofa(covm,ldcovm,p,info) IF (.NOT. (info.NE.0)) GO TO 30 WRITE (*,*) ' COVM not positive definite in SETGMN' - CALL XSTOPX (' COVM not positive definite in SETGMN') + STOP ' COVM not positive definite in SETGMN' 30 icount = p + 1 C
--- a/libcruft/ranlib/setsd.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/setsd.f Thu Oct 15 06:02:21 1998 +0000 @@ -62,8 +62,7 @@ IF (qrgnin()) GO TO 10 WRITE (*,*) ' SETSD called before random number generator ', + ' initialized -- abort!' - CALL XSTOPX - + (' SETSD called before random number generator initialized') + STOP ' SETSD called before random number generator initialized' 10 CALL getcgn(g) ig1(g) = iseed1
--- a/libcruft/ranlib/sexpo.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/sexpo.f Thu Oct 15 06:02:21 1998 +0000 @@ -23,14 +23,32 @@ C C C**********************************************************************C C - DIMENSION q(8) - EQUIVALENCE (q(1),q1) C C Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N C (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION C +C JJV added a Save statement for q (in Data statement) +C .. Local Scalars .. + REAL a,q1,u,umin,ustar + INTEGER i +C .. +C .. Local Arrays .. + REAL q(8) +C .. +C .. External Functions .. + REAL ranf + EXTERNAL ranf +C .. +C .. Equivalences .. + EQUIVALENCE (q(1),q1) +C .. +C .. Save statement .. + SAVE q +C .. +C .. Data statements .. DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833, + .9999986,.9999999/ +C .. C 10 a = 0.0 u = ranf() @@ -38,7 +56,10 @@ 20 a = a + q1 30 u = u + u - IF (u.LE.1.0) GO TO 20 +C JJV changed the following to reflect the true algorithm and +C JJV prevent unpredictable behavior if U is initially 0.5. +C IF (u.LE.1.0) GO TO 20 + IF (u.LT.1.0) GO TO 20 40 u = u - 1.0 IF (u.GT.q1) GO TO 60 50 sexpo = a + u
--- a/libcruft/ranlib/sgamma.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/sgamma.f Thu Oct 15 06:02:21 1998 +0000 @@ -50,19 +50,41 @@ C COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K) C COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K) C +C .. Scalar Arguments .. + REAL a +C .. +C .. Local Scalars .. (JJV added B0 to fix rare and subtle bug) + REAL a1,a2,a3,a4,a5,a6,a7,aa,aaa,b,b0,c,d,e,e1,e2,e3,e4,e5,p,q,q0, + + q1,q2,q3,q4,q5,q6,q7,r,s,s2,si,sqrt32,t,u,v,w,x +C .. +C .. External Functions .. + REAL ranf,sexpo,snorm + EXTERNAL ranf,sexpo,snorm +C .. +C .. Intrinsic Functions .. + INTRINSIC abs,alog,exp,sign,sqrt +C .. +C .. Save statement .. +C JJV added Save statement for vars in Data satatements + SAVE aa,aaa,s2,s,d,q0,b,si,c,q1,q2,q3,q4,q5,q6,q7,a1,a2,a3,a4,a5, + + a6,a7,e1,e2,e3,e4,e5,sqrt32 +C .. +C .. Data statements .. +C +C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" +C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 +C DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121, + -.00007388,.00024511,.00024240/ DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921, + .1423657,-.1367177,.1233795/ DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/ -C -C PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A" -C SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380 -C DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/ +C .. +C .. Executable Statements .. C IF (a.EQ.aa) GO TO 10 - IF (a.LT.1.0) GO TO 120 + IF (a.LT.1.0) GO TO 130 C C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED C @@ -115,9 +137,9 @@ C C CASE 1: A .LE. 3.686 C - 30 b = .463 + s - .178*s2 + 30 b = .463 + s + .178*s2 si = 1.235 - c = .195/s - .079 + .016*s + c = .195/s - .079 + .16*s C C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE C @@ -147,45 +169,67 @@ C C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 C - IF (t.LT. (-.7187449)) GO TO 70 + 80 IF (t.LT. (-.7187449)) GO TO 70 C C STEP 10: CALCULATION OF V AND QUOTIENT Q C v = t/ (s+s) - IF (abs(v).LE.0.25) GO TO 80 + IF (abs(v).LE.0.25) GO TO 90 q = q0 - s*t + 0.25*t*t + (s2+s2)*alog(1.0+v) - GO TO 90 + GO TO 100 - 80 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v + 90 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v C C STEP 11: HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8) C - 90 IF (q.LE.0.0) GO TO 70 - IF (q.LE.0.5) GO TO 100 - w = exp(q) - 1.0 - GO TO 110 + 100 IF (q.LE.0.0) GO TO 70 + IF (q.LE.0.5) GO TO 110 +C +C JJV modified the code through line 125 to handle large Q case +C + IF (q.LT.15.0) GO TO 105 +C +C JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for real Q) +C JJV so reformulate test at 120 in terms of one EXP, if not too big +C JJV 87.49823 is close to the largest real which can be +C JJV exponentiated (87.49823 = log(1.0E38)) +C + IF ((q+e-0.5*t*t).GT.87.49823) GO TO 125 + IF (c*abs(u).GT.exp(q+e-0.5*t*t)) GO TO 70 + GO TO 125 - 100 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q + 105 w = exp(q) - 1.0 + GO TO 120 + + 110 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q C C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 C - 110 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70 - x = s + 0.5*t + 120 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70 + 125 x = s + 0.5*t sgamma = x*x RETURN C C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) C - 120 aa = 0.0 - b = 1.0 + .3678794*a - 130 p = b*ranf() - IF (p.GE.1.0) GO TO 140 +C JJV changed B to B0 (which was added to declarations for this) +C JJV in 130 to END to fix rare and subtle bug. +C JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful). +C JJV Reasons: the state of AA only serves to tell the A .GE. 1.0 +C JJV case if certain A-dependant constants need to be recalculated. +C JJV The A .LT. 1.0 case (here) no longer changes any of these, and +C JJV the recalculation of B (which used to change with an +C JJV A .LT. 1.0 call) is governed by the state of AAA anyway. +C + 130 b0 = 1.0 + .3678794*a + 140 p = b0*ranf() + IF (p.GE.1.0) GO TO 150 sgamma = exp(alog(p)/a) - IF (sexpo().LT.sgamma) GO TO 130 + IF (sexpo().LT.sgamma) GO TO 140 RETURN - 140 sgamma = -alog((b-p)/a) - IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 130 + 150 sgamma = -alog((b0-p)/a) + IF (sexpo().LT. (1.0-a)*alog(sgamma)) GO TO 140 RETURN END
--- a/libcruft/ranlib/snorm.f Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/snorm.f Thu Oct 15 06:02:21 1998 +0000 @@ -23,11 +23,29 @@ C C C**********************************************************************C C - DIMENSION a(32),d(31),t(31),h(31) C C THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND C H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE C +C .. Local Scalars .. + REAL aa,s,tt,u,ustar,w,y + INTEGER i +C .. +C .. Local Arrays .. + REAL a(32),d(31),h(31),t(31) +C .. +C .. External Functions .. + REAL ranf + EXTERNAL ranf +C .. +C .. Intrinsic Functions .. + INTRINSIC float,int +C .. +C .. Save statement .. +C JJV added a Save statement for arrays initialized in Data statmts + SAVE a,d,t,h +C .. +C .. Data statements .. DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991, + .2372021,.2776904,.3186394,.3601299,.4022501,.4450965, + .4887764,.5334097,.5791322,.6260990,.6744898,.7245144, @@ -53,6 +71,8 @@ + .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1, + .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016, + .7010474/ +C .. +C .. Executable Statements .. C 10 u = ranf() s = 0.0
--- a/libcruft/ranlib/tstgmn.for Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/tstgmn.for Thu Oct 15 06:02:21 1998 +0000 @@ -1,4 +1,7 @@ - REAL FUNCTION covar(x,y,n) +C JJV changed name to ONECOV to avoid confusion with array COVAR +C JJV this was also changed in the body of the function +C REAL FUNCTION covar(x,y,n) + REAL FUNCTION onecov(x,y,n) C .. Scalar Arguments .. INTEGER n C .. @@ -18,21 +21,33 @@ C .. Executable Statements .. CALL stat(x,n,avx,varx,xmin,xmax) CALL stat(y,n,avy,vary,xmin,xmax) - covar = 0.0 +C covar = 0.0 + onecov = 0.0 DO 10,i = 1,n - covar = covar + (x(i)-avx)* (y(i)-avy) - 10 CONTINUE - covar = covar/real(n-1) +C covar = covar + (x(i)-avx)* (y(i)-avy) + onecov = onecov + (x(i)-avx)* (y(i)-avy) + 10 CONTINUE +C covar = covar/real(n-1) + onecov = onecov/real(n-1) RETURN - + END - SUBROUTINE prcomp(p,mean,xcovar,answer) - INTEGER p,maxp +C JJV Added argument LDXCOV (leading dimension of XCOVAR) to be +C JJV consistent with the program TSTGMN, see comments below. +C JJV This change necessitated changes in the declarations. +C SUBROUTINE prcomp(p,mean,xcovar,answer) + SUBROUTINE prcomp(p,mean,xcovar,ldxcov,answer) + +C INTEGER p,maxp + INTEGER p,maxp,ldxcov PARAMETER (maxp=10) - REAL mean(p),xcovar(p,p),rcovar(maxp,maxp) +C REAL mean(p),xcovar(p,p),rcovar(maxp,maxp) + REAL mean(p),xcovar(ldxcov,p),rcovar(maxp,maxp) REAL answer(1000,maxp) - REAL rmean(maxp),rvar(maxp) +C JJV added ONECOV because of name change to function COVAR +C REAL rmean(maxp),rvar(maxp) + REAL rmean(maxp),rvar(maxp),onecov INTEGER maxobs PARAMETER (maxobs=1000) @@ -46,7 +61,9 @@ DO 30,i = 1,p DO 20,j = 1,i - 1 WRITE (*,*) ' I = ',i,' J = ',j - rcovar(i,j) = covar(answer(1,i),answer(1,j),maxobs) +C JJV changed COVAR to match new name +C rcovar(i,j) = covar(answer(1,i),answer(1,j),maxobs) + rcovar(i,j) = onecov(answer(1,i),answer(1,j),maxobs) WRITE (*,*) ' Covariance ',xcovar(i,j),' Generated ', + rcovar(i,j) 20 CONTINUE @@ -54,14 +71,21 @@ RETURN END - SUBROUTINE setcov(p,var,corr,covar) + +C JJV added LDCOV (leading dimension of COVAR) to be +C JJV consistent with the program TSTGMN, see comments below. +C JJV This change necessitated changes in the declarations. +C SUBROUTINE setcov(p,var,corr,covar) + SUBROUTINE setcov(p,var,corr,covar,ldcov) C Set covariance matrix from variance and common correlation C .. Scalar Arguments .. REAL corr - INTEGER p +C INTEGER p + INTEGER p,ldcov C .. C .. Array Arguments .. - REAL covar(p,p),var(p) +C REAL covar(p,p),var(p) + REAL covar(ldcov,p),var(p) C .. C .. Local Scalars .. INTEGER i,j @@ -83,6 +107,7 @@ RETURN END + SUBROUTINE stat(x,n,av,var,xmin,xmax) C .. Scalar Arguments .. REAL av,var,xmax,xmin @@ -116,15 +141,24 @@ RETURN END + PROGRAM tstgmn C Test Generation of Multivariate Normal Data +C JJV SETGMN was: SUBROUTINE setgmn(meanv,covm,p,parm) +C JJV is: SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm) +C JJV So the covariance matrices have been changed to 2-dim'l +C JJV matrices, and the additional argument has been added to +C JJV the subroutine call. Additional changes have been made +C JJV to reflect this. (in declarations, the matrix copy routine, +C JJV and in subroutine calls.) C .. Parameters .. INTEGER maxp PARAMETER (maxp=10) INTEGER maxobs PARAMETER (maxobs=1000) - INTEGER p2 - PARAMETER (p2=maxp*maxp) +C JJV this parameter is no longer needed +C INTEGER p2 +C PARAMETER (p2=maxp*maxp) C .. C .. Local Scalars .. REAL corr @@ -132,8 +166,10 @@ CHARACTER phrase*100 C .. C .. Local Arrays .. - REAL answer(1000,maxp),ccovar(p2),covar(p2),mean(maxp),param(500), - + temp(maxp),var(maxp),work(maxp) +C REAL answer(1000,maxp),ccovar(p2),covar(p2),mean(maxp),param(500), +C + temp(maxp),var(maxp),work(maxp) + REAL answer(1000,maxp),ccovar(maxp,maxp),covar(maxp,maxp), + + mean(maxp),param(500),temp(maxp),var(maxp),work(maxp) C .. C .. External Subroutines .. EXTERNAL genmn,phrtsd,prcomp,setall,setcov,setgmn @@ -158,25 +194,33 @@ READ (*,*) (var(i),i=1,p) WRITE (*,*) 'Enter correlation of all variables' READ (*,*) corr - CALL setcov(p,var,corr,covar) +C CALL setcov(p,var,corr,covar) + CALL setcov(p,var,corr,covar,maxp) WRITE (*,*) ' Enter phrase to initialize rn generator' READ (*,'(a)') phrase CALL phrtsd(phrase,is1,is2) CALL setall(is1,is2) - DO 20,i = 1,p2 - ccovar(i) = covar(i) - 20 CONTINUE +C DO 20,i = 1,p2 +C ccovar(i) = covar(i) +C 20 CONTINUE + DO 25,i = 1,maxp + DO 20,j = 1,maxp + ccovar(i,j) = covar(i,j) + 20 CONTINUE + 25 CONTINUE C C Generate Variables C - CALL setgmn(mean,ccovar,p,param) +C CALL setgmn(mean,ccovar,p,param) + CALL setgmn(mean,ccovar,maxp,p,param) DO 40,iobs = 1,maxobs CALL genmn(param,work,temp) DO 30,j = 1,p answer(iobs,j) = work(j) 30 CONTINUE 40 CONTINUE - CALL prcomp(p,mean,covar,answer) +C CALL prcomp(p,mean,covar,answer) + CALL prcomp(p,mean,covar,maxp,answer) C C Print Comparison of Generated and Reconstructed Values C
--- a/libcruft/ranlib/tstmid.for Fri Oct 09 03:24:56 1998 +0000 +++ b/libcruft/ranlib/tstmid.for Thu Oct 15 06:02:21 1998 +0000 @@ -48,28 +48,25 @@ IMPLICIT LOGICAL (q) C Interactive test for PHRTSD C .. Parameters .. - INTEGER mxwh - PARAMETER (mxwh=10) + INTEGER mxwh,mxncat + PARAMETER (mxwh=15,mxncat=100) C .. C .. Local Scalars .. - REAL av,avtr,var,vartr,xmin,xmax - INTEGER i,is1,is2,itmp,iwhich,j,mxint,nperm,nrep,ntot + REAL av,avtr,var,vartr,xmin,xmax,pevt,psum,rtry + INTEGER i,is1,is2,itmp,iwhich,j,mxint,nperm,nrep,ntot,ntry,ncat CHARACTER type*4,phrase*100 C .. C .. Local Arrays .. - REAL array(1000),param(3) + REAL array(1000),param(3),prob(mxncat) INTEGER iarray(1000),perm(500) C .. C .. External Functions .. - REAL genbet,genchi,genf,gennch,gennf,genunf - INTEGER ignuin,myhand - EXTERNAL genbet,genchi,genf,gennch,gennf,genunf,ignuin,myhand + REAL genbet,genchi,genf,gennch,gennf,genunf,genexp,gengam,gennor + INTEGER ignuin,ignnbn + EXTERNAL genbet,genchi,genf,gennch,gennf,genunf,ignuin C .. C .. External Subroutines .. - EXTERNAL genprm,phrtsd,setall,stat,trstat -C .. -C .. Equivalences .. - EQUIVALENCE (array,iarray) + EXTERNAL genprm,phrtsd,setall,stat,trstat,genmul C .. C .. Executable Statements .. WRITE (*,9000) @@ -82,7 +79,15 @@ + ' and prints it.'/ + ' For uniform integers asks for upper bound, number of'/ + ' replicates per integer in 1..upper bound.'/ - + ' Prints table of num times each integer generated.') + + ' Prints table of num times each integer generated.'/ + + ' For multinomial asks for number of events to be'/ + + ' classified, number of categories in which they'/ + + ' are to be classified, and the probabilities that'/ + + ' an event will be classified in the categories,'/ + + ' for all but the last category. Prints table of'/ + + ' number of events by category, true probability'/ + + ' associated with each category, and observed'/ + + ' proportion of events in each category.') C C Menu for choosing tests C @@ -99,7 +104,12 @@ + ' (7) Generate uniform reals'/ + ' (8) Generate beta deviates'/ + ' (9) Generate binomial outcomes'/ - + ' (10) Generate Poisson ourcomes'/) + + ' (10) Generate Poisson outcomes'/ + + ' (11) Generate exponential deviates'/ + + ' (12) Generate gamma deviates'/ + + ' (13) Generate multinomial outcomes'/ + + ' (14) Generate normal deviates'/ + + ' (15) Generate negative binomial outcomes'/) READ (*,*) iwhich IF (.NOT. (iwhich.LT.0.OR.iwhich.GT.mxwh)) GO TO 20 @@ -128,9 +138,9 @@ 9020 FORMAT (' Mean Generated: ',T30,G15.7,5X,'True:',T60, + G15.7/' Variance Generated:',T30,G15.7,5X,'True:',T60, - + G15.7/' Minimun: ',T30,G15.7,5X,'Maximum:',T60,G15.7) + + G15.7/' Minimum: ',T30,G15.7,5X,'Maximum:',T60,G15.7) - GO TO 280 + GO TO 420 40 IF ((2).NE. (iwhich)) GO TO 60 @@ -147,7 +157,7 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 60 IF ((3).NE. (iwhich)) GO TO 80 @@ -164,11 +174,10 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 80 IF ((4).NE. (iwhich)) GO TO 100 - C C Noncentral F deviates C @@ -183,7 +192,7 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 100 IF ((5).NE. (iwhich)) GO TO 140 @@ -204,7 +213,7 @@ CALL genprm(perm,nperm) WRITE (*,*) ' Perm Generated' WRITE (*,'(20I4)') (perm(i),i=1,nperm) - GO TO 280 + GO TO 420 140 IF ((6).NE. (iwhich)) GO TO 170 @@ -225,7 +234,7 @@ 160 CONTINUE WRITE (*,*) ' Counts of Integers Generated' WRITE (*,'(20I4)') (iarray(j),j=1,mxint) - GO TO 280 + GO TO 420 170 IF ((7).NE. (iwhich)) GO TO 190 @@ -241,7 +250,7 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 190 IF ((8).NE. (iwhich)) GO TO 210 @@ -257,9 +266,10 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 210 IF ((9).NE. (iwhich)) GO TO 240 + C C Binomial outcomes C @@ -278,9 +288,10 @@ param(2) = pevt CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 240 IF ((10).NE. (iwhich)) GO TO 270 + C C Poisson outcomes C @@ -296,10 +307,126 @@ CALL stat(array,1000,av,var,xmin,xmax) CALL trstat(type,param,avtr,vartr) WRITE (*,9020) av,avtr,var,vartr,xmin,xmax - GO TO 280 + GO TO 420 + + 270 IF ((11).NE. (iwhich)) GO TO 290 + +C +C Exponential deviates +C + type = 'expo' + WRITE (*,*) ' Enter (real) AV for Exponential' + READ (*,*) param(1) + DO 280,i = 1,1000 + array(i) = genexp(param(1)) + 280 CONTINUE + CALL stat(array,1000,av,var,xmin,xmax) + CALL trstat(type,param,avtr,vartr) + WRITE (*,9020) av,avtr,var,vartr,xmin,xmax + + GO TO 420 + + 290 IF ((12).NE. (iwhich)) GO TO 310 + +C +C Gamma deviates +C + type = 'gamm' + WRITE (*,*) ' Enter (real) A, (real) R for Gamma deviate' + READ (*,*) param(1),param(2) + DO 300,i = 1,1000 + array(i) = gengam(param(1),param(2)) + 300 CONTINUE + CALL stat(array,1000,av,var,xmin,xmax) + CALL trstat(type,param,avtr,vartr) + WRITE (*,9020) av,avtr,var,vartr,xmin,xmax + GO TO 420 + + 310 IF ((13).NE. (iwhich)) GO TO 360 - 270 CONTINUE - 280 GO TO 10 +C +C Multinomial outcomes +C + WRITE (*,*) ' Enter (int) number of observations: ' + READ (*,*) ntry + 320 WRITE (*,*) ' Enter (int) num. of categories: <= ',mxncat + READ (*,*) ncat + IF (ncat.GT.mxncat) THEN + WRITE (*,*) ' number of categories must be <= ',mxncat + WRITE (*,*) ' Try again ... ' + GO TO 320 + END IF + WRITE (*,*) ' Enter (real) prob. vector of length ',ncat-1 + READ (*,*) (prob(i),i=1,ncat-1) + CALL genmul(ntry,prob,ncat,iarray) + ntot = 0 + IF (ntry.GT.0) THEN + rtry = real(ntry) + DO 330, i = 1,ncat + ntot = ntot + iarray(i) + array(i) = iarray(i)/rtry + 330 CONTINUE + ELSE + DO 340, i = 1,ncat + ntot = ntot + iarray(i) + array(i) = 0.0 + 340 CONTINUE + ENDIF + psum = 0.0 + DO 350, i = 1,ncat-1 + psum = psum + prob(i) + 350 CONTINUE + prob(ncat) = 1.0 - psum + + WRITE (*,*) ' Total number of observations: ',ntot + WRITE (*,*) ' Total observations by category: ' + WRITE (*,'(10I8)') (iarray(i),i=1,ncat) + WRITE (*,*) ' True probabilities by category: ' + WRITE (*,'(8F10.7)') (prob(i),i=1,ncat) + WRITE (*,*) ' Observed proportions by category: ' + WRITE (*,'(8F10.7)') (array(i),i=1,ncat) + GO TO 420 + + 360 IF ((14).NE. (iwhich)) GO TO 380 + +C +C Normal deviates +C + type = 'norm' + WRITE (*,*) ' Enter (real) AV, (real) SD for Normal' + READ (*,*) param(1),param(2) + DO 370,i = 1,1000 + array(i) = gennor(param(1),param(2)) + 370 CONTINUE + CALL stat(array,1000,av,var,xmin,xmax) + CALL trstat(type,param,avtr,vartr) + WRITE (*,9020) av,avtr,var,vartr,xmin,xmax + GO TO 420 + + 380 IF ((15).NE. (iwhich)) GO TO 410 + +C +C Negative Binomial outcomes +C + type = 'nbin' + WRITE (*,*) ' Enter required (int) Number of events then ' + WRITE (*,*) ' (real) Prob of an event for negative binomial' + READ (*,*) ntry,pevt + DO 390,i = 1,1000 + iarray(i) = ignnbn(ntry,pevt) + 390 CONTINUE + DO 400,i = 1,1000 + array(i) = iarray(i) + 400 CONTINUE + CALL stat(array,1000,av,var,xmin,xmax) + param(1) = ntry + param(2) = pevt + CALL trstat(type,param,avtr,vartr) + WRITE (*,9020) av,avtr,var,vartr,xmin,xmax + GO TO 420 + + 410 CONTINUE + 420 GO TO 10 END SUBROUTINE trstat(type,parin,av,var) @@ -323,6 +450,12 @@ C 'ncf' noncentral f C 'unif' uniform C 'beta' beta distribution +C 'bin' binomial +C 'pois' poisson +C 'expo' exponential +C 'gamm' gamma +C 'norm' normal +C 'nbin' negative binomial C CHARACTER*(4) TYPE C C PARIN --> Array containing parameters of distribution @@ -344,12 +477,23 @@ C beta C PARIN(1) is A C PARIN(2) is B -C REAL PARIN(*) -C binonial +C binomial C PARIN(1) is Number of trials C PARIN(2) is Prob Event at Each Trial C poisson C PARIN(1) is Mean +C exponential +C PARIN(1) is Mean +C gamma +C PARIN(1) is A +C PARIN(2) is R +C normal +C PARIN(1) is Mean +C PARIN(2) is Standard Deviation +C negative binomial +C PARIN(1) is required Number of events +C PARIN(2) is Probability of event +C REAL PARIN(*) C C AV <-- Mean of specified distribution with specified parameters C REAL AV @@ -378,14 +522,14 @@ IF (('chis').NE. (type)) GO TO 10 av = parin(1) var = 2.0*parin(1) - GO TO 170 + GO TO 210 10 IF (('ncch').NE. (type)) GO TO 20 a = parin(1) + parin(2) b = parin(2)/a av = a var = 2.0*a* (1.0+b) - GO TO 170 + GO TO 210 20 IF (('f').NE. (type)) GO TO 70 IF (.NOT. (parin(2).LE.2.0001)) GO TO 30 @@ -399,7 +543,7 @@ 50 var = (2.0*parin(2)**2* (parin(1)+parin(2)-2.0))/ + (parin(1)* (parin(2)-2.0)**2* (parin(2)-4.0)) - 60 GO TO 170 + 60 GO TO 210 70 IF (('ncf').NE. (type)) GO TO 120 IF (.NOT. (parin(2).LE.2.0001)) GO TO 80 @@ -415,34 +559,53 @@ + (parin(2)-2.0) b = (parin(2)-2.0)**2* (parin(2)-4.0) var = 2.0* (parin(2)/parin(1))**2* (a/b) - 110 GO TO 170 + 110 GO TO 210 120 IF (('unif').NE. (type)) GO TO 130 range = parin(2) - parin(1) av = parin(1) + range/2.0 var = range**2/12.0 - GO TO 170 + GO TO 210 130 IF (('beta').NE. (type)) GO TO 140 av = parin(1)/ (parin(1)+parin(2)) var = (av*parin(2))/ ((parin(1)+parin(2))* + (parin(1)+parin(2)+1.0)) - WRITE (*,*) ' A, B, AV, VAR ',parin(1),parin(2),av,var - GO TO 170 + GO TO 210 140 IF (('bin').NE. (type)) GO TO 150 av = parin(1)*parin(2) var = av* (1.0-parin(2)) - GO TO 170 + GO TO 210 150 IF (('pois').NE. (type)) GO TO 160 av = parin(1) var = parin(1) - GO TO 170 + GO TO 210 + + 160 IF (('expo').NE. (type)) GO TO 170 + av = parin(1) + var = parin(1)**2 + GO TO 210 + + 170 IF (('gamm').NE. (type)) GO TO 180 + av = parin(2) / parin(1) + var = av / parin(1) + GO TO 210 - 160 WRITE (*,*) 'Unimplemented type ',type + 180 IF (('norm').NE. (type)) GO TO 190 + av = parin(1) + var = parin(2)**2 + GO TO 210 + + 190 IF (('nbin').NE. (type)) GO TO 200 + av = parin(1) * (1.0 - parin(2)) / parin(2) + var = av / parin(2) + GO TO 210 + + 200 WRITE (*,*) 'Unimplemented type ',type STOP 'Unimplemented type in TRSTAT' - 170 RETURN + 210 RETURN END