SUBROUTINE ADVNST(K) 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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. INTEGER K C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G,I,IB1,IB2 C .. C .. EXTERNAL FUNCTIONS .. INTEGER MLTMOD LOGICAL QRGNIN EXTERNAL MLTMOD,QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN,SETSD C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C ABORT UNLESS RANDOM NUMBER GENERATOR INITIALIZED IF (QRGNIN()) GO TO 10 WRITE (*,*) ' ADVNST CALLED BEFORE RANDOM NUMBER GENERATOR ', + ' INITIALIZED -- ABORT!' STOP ' ADVNST CALLED BEFORE RANDOM NUMBER GENERATOR INITIALIZED' 10 CALL GETCGN(G) C IB1 = A1 IB2 = A2 DO 20,I = 1,K IB1 = MLTMOD(IB1,IB1,M1) IB2 = MLTMOD(IB2,IB2,M2) 20 CONTINUE CALL SETSD(MLTMOD(IB1,CG1(G),M1),MLTMOD(IB2,CG2(G),M2)) C C NOW, IB1 = A1**K AND IB2 = A2**K C RETURN END C------------------------------------------------------------------ DOUBLE PRECISION FUNCTION GENBET(AA,BB) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION A C C B --> SECOND PARAMETER OF THE BETA DISTRIBUTION C DOUBLE PRECISION B C C C METHOD C C C R. C. H. CHENG C GENERATING BETA VARIATEW WITH NONINTEGRAL SHAPE PARAMETERS C COMMUNICATIONS OF THE ACM, 21:317-322 (1978) C (ALGORITHMS BB AND BC) C C********************************************************************** C .. PARAMETERS .. C CLOSE TO THE LARGEST NUMBER THAT CAN BE EXPONENTIATED DOUBLE PRECISION EXPMAX PARAMETER (EXPMAX=89.0D0) C CLOSE TO THE LARGEST REPRESENTABLE SINGLE PRECISION NUMBER DOUBLE PRECISION INFNTY PARAMETER (INFNTY=1.0D38) C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION AA,BB C .. C .. LOCAL SCALARS .. DOUBLE PRECISION A,ALPHA,B,BETA,DELTA,GAMMA,K1,K2,OLDA,OLDB,R, +S,T,U1,U2,V,W,Y, Z LOGICAL QSAME C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION RANF EXTERNAL RANF C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DEXP,DLOG,DMAX1,DMIN1,DSQRT C .. C .. SAVE STATEMENT .. SAVE OLDA,OLDB,ALPHA,BETA,GAMMA,K1,K2 C .. C .. DATA STATEMENTS .. DATA OLDA,OLDB/-1,-1/ 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!' WRITE (*,*) ' AA: ',AA,' BB ',BB STOP ' AA OR BB <= 0 IN GENBET - ABORT!' 10 OLDA = AA OLDB = BB 20 IF (.NOT. (DMIN1(AA,BB).GT.1.0D0)) GO TO 100 C ALBORITHM BB C C INITIALIZE C IF (QSAME) GO TO 30 A = MIN(AA,BB) B = MAX(AA,BB) ALPHA = A + B BETA = DSQRT((ALPHA-2.0D0)/ (2.0D0*A*B-ALPHA)) GAMMA = A + 1.0/BETA 30 CONTINUE 40 U1 = RANF() C C STEP 1 C U2 = RANF() V = BETA*DLOG(U1/ (1.0D0-U1)) IF (.NOT. (V.GT.EXPMAX)) GO TO 50 W = INFNTY GO TO 60 50 W = A*DEXP(V) 60 Z = U1**2D0*U2 R = GAMMA*V - 1.3862944D0 S = A + R - W C C STEP 2 C IF ((S+2.609438D0).GE. (5.0D0*Z)) GO TO 70 C C STEP 3 C T = DLOG(Z) IF (S.GT.T) GO TO 70 C C STEP 4 C IF ((R+ALPHA*DLOG(ALPHA/ (B+W))).LT.T) GO TO 40 C C STEP 5 C 70 IF (.NOT. (AA.EQ.A)) GO TO 80 GENBET = W/ (B+W) GO TO 90 80 GENBET = B/ (B+W) 90 GO TO 230 C ALGORITHM BC C C INITIALIZE C 100 IF (QSAME) GO TO 110 A =DMAX1(AA,BB) B = DMIN1(AA,BB) ALPHA = A + B BETA = 1.0/B DELTA = 1.0D0 + A - B K1 = DELTA* (0.0138889D0+0.0416667D0*B)/ (A*BETA-0.777778D0) K2 = 0.25D0 + (0.5D0+0.25D0/DELTA)*B 110 CONTINUE 120 U1 = RANF() C C STEP 1 C U2 = RANF() IF (U1.GE.0.5) GO TO 130 C C STEP 2 C Y = U1*U2 Z = U1*Y IF ((0.25D0*U2+Z-Y).GE.K1) GO TO 120 GO TO 170 C C STEP 3 C 130 Z = U1**2.0D0*U2 IF (.NOT. (Z.LE.0.25D0)) GO TO 160 V = BETA*DLOG(U1/ (1.0D0-U1)) IF (.NOT. (V.GT.EXPMAX)) GO TO 140 W = INFNTY GO TO 150 140 W = A*DEXP(V) 150 GO TO 200 160 IF (Z.GE.K2) GO TO 120 C C STEP 4 C C C STEP 5 C 170 V = BETA*DLOG(U1/ (1.0D0-U1)) IF (.NOT. (V.GT.EXPMAX)) GO TO 180 W = INFNTY GO TO 190 180 W = A*DEXP(V) 190 IF ((ALPHA* (DLOG(ALPHA/ (B+W))+V)-1.3862944D0).LT.DLOG(Z)) +GO TO 120 C C STEP 6 C 200 IF (.NOT. (A.EQ.AA)) GO TO 210 GENBET = W/ (B+W) GO TO 220 210 GENBET = B/ (B+W) 220 CONTINUE 230 RETURN END C--------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENCHI(DF) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION DF C C C METHOD C C C USES RELATION BETWEEN CHISQUARE AND GAMMA. C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GENGAM EXTERNAL GENGAM C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (DF.LE.0.0D0)) GO TO 10 WRITE (*,*) 'DF <= 0 IN GENCHI - ABORT' WRITE (*,*) 'VALUE OF DF: ',DF STOP 'DF <= 0 IN GENCHI - ABORT' 10 GENCHI = 2.0D0*GENGAM(1.0D0,DF/2.0D0) RETURN END C---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENEXP(AV) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION AV C C GENEXP <-- THE RANDOM DEVIATE. C DOUBLE PRECISION 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 .. SCALAR ARGUMENTS .. DOUBLE PRECISION AV C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION SEXPO EXTERNAL SEXPO C .. C .. EXECUTABLE STATEMENTS .. GENEXP = SEXPO()*AV RETURN END C-------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENF(DFN,DFD) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION DFN C DFD --> DENOMINATOR DEGREES OF FREEDOM C (MUST BE POSITIVE) C DOUBLE PRECISION DFD C C C METHOD C C C DIRECTLY GENERATES RATIO OF CHISQUARE VARIATES C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION DFD,DFN C .. C .. LOCAL SCALARS .. DOUBLE PRECISION XDEN,XNUM C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GENCHI EXTERNAL GENCHI C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (DFN.LE.0.0D0.OR.DFD.LE.0.0D0)) GO TO 10 WRITE (*,*) 'DEGREES OF FREEDOM NONPOSITIVE IN GENF - ABORT!' WRITE (*,*) 'DFN VALUE: ',DFN,'DFD VALUE: ',DFD STOP 'DEGREES OF FREEDOM NONPOSITIVE IN GENF - ABORT!' 10 XNUM = GENCHI(DFN)/DFN C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) XDEN = GENCHI(DFD)/DFD IF (.NOT. (XDEN.LE. (1.0D-38*XNUM))) GO TO 20 WRITE (*,*) ' GENF - GENERATED NUMBERS WOULD CAUSE OVERFLOW' WRITE (*,*) ' NUMERATOR ',XNUM,' DENOMINATOR ',XDEN WRITE (*,*) ' GENF RETURNING 1.0E38' GENF = 1.0D38 GO TO 30 20 GENF = XNUM/XDEN 30 RETURN END C-------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENGAM(A,R) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION A C C R --> SHAPE PARAMETER OF GAMMA DISTRIBUTION C DOUBLE PRECISION 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 .. SCALAR ARGUMENTS .. DOUBLE PRECISION A,R C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION SGAMMA EXTERNAL SGAMMA C .. C .. EXECUTABLE STATEMENTS .. GENGAM = SGAMMA(R) GENGAM = GENGAM/A RETURN END C--------------------------------------------------------------------- SUBROUTINE GENMN(PARM,X,WORK) 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 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 DOUBLE PRECISION PARM(*) C C X <-- VECTOR DEVIATE GENERATED. C DOUBLE PRECISION X(P) C C WORK <--> SCRATCH ARRAY C DOUBLE PRECISION WORK(P) C C C METHOD C C C 1) GENERATE P INDEPENDENT STANDARD NORMAL DEVIATES - EI ~ N(0,1) C C 2) USING CHOLESKY DECOMPOSITION FIND A S.T. TRANS(A)*A = COVM C C 3) TRANS(A)E + MEANV ~ N(MEANV,COVM) C C********************************************************************** C .. ARRAY ARGUMENTS .. DOUBLE PRECISION PARM(*),WORK(*),X(*) C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AE INTEGER I,ICOUNT,J,P C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION SNORM EXTERNAL SNORM C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC IDINT C .. C .. EXECUTABLE STATEMENTS .. P = IDINT(PARM(1)) C C GENERATE P INDEPENDENT NORMAL DEVIATES - WORK ~ N(0,1) C DO 10,I = 1,P WORK(I) = SNORM() 10 CONTINUE DO 30,I = 1,P C C PARM (P+2 : P*(P+3)/2 + 1) CONTAINS A, THE CHOLESKY C DECOMPOSITION OF THE DESIRED COVARIANCE MATRIX. C TRANS(A)(1,1) = PARM(P+2) C TRANS(A)(2,1) = PARM(P+3) C TRANS(A)(2,2) = PARM(P+2+P) C TRANS(A)(3,1) = PARM(P+4) C TRANS(A)(3,2) = PARM(P+3+P) C TRANS(A)(3,3) = PARM(P+2-1+2P) ... C C TRANS(A)*WORK + MEANV ~ N(MEANV,COVM) C ICOUNT = 0 AE = 0.0D0 DO 20,J = 1,I ICOUNT = ICOUNT + J - 1 AE = AE + PARM(I+ (J-1)*P-ICOUNT+P+1)*WORK(J) 20 CONTINUE X(I) = AE + PARM(I+1) 30 CONTINUE RETURN C END C--------------------------------------------------------------------- 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 DOUBLE PRECISION 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 .. DOUBLE PRECISION P(*) INTEGER IX(*) C .. C .. LOCAL SCALARS .. DOUBLE PRECISION PROB,PTOT,SUM INTEGER I,ICAT,NTOT C .. C .. EXTERNAL FUNCTIONS .. INTEGER IGNBIN EXTERNAL IGNBIN C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS 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.0D0) STOP 'SOME P(I) < 0 IN GENMUL' IF (P(I).GT.1.0D0) 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 C------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION GENNCH(DF,XNONC) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION DF C C XNONC --> NONCENTRALITY PARAMETER OF THE CHISQUARE C (MUST BE >= 0.0) C DOUBLE PRECISION 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 .. SCALAR ARGUMENTS .. DOUBLE PRECISION DF,XNONC C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GENCHI,GENNOR EXTERNAL GENCHI,GENNOR C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DSQRT C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (DF.LE.1.0D0.OR.XNONC.LT.0.0D0)) GO TO 10 WRITE (*,*) 'DF <= 1 OR XNONC < 0 IN GENNCH - ABORT' WRITE (*,*) 'VALUE OF DF: ',DF,' VALUE OF XNONC',XNONC STOP 'DF <= 1 OR XNONC < 0 IN GENNCH - ABORT' 10 GENNCH = GENCHI(DF-1.0D0) + GENNOR(DSQRT(XNONC),1.0D0)**2 RETURN END C---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENNF(DFN,DFD,XNONC) C********************************************************************** C C DOUBLE PRECISION 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 DOUBLE PRECISION DFN C DFD --> DENOMINATOR DEGREES OF FREEDOM C (MUST BE POSITIVE) C DOUBLE PRECISION DFD C C XNONC --> NONCENTRALITY PARAMETER C (MUST BE NONNEGATIVE) C DOUBLE PRECISION 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 .. SCALAR ARGUMENTS .. DOUBLE PRECISION DFD,DFN,XNONC C .. C .. LOCAL SCALARS .. DOUBLE PRECISION XDEN,XNUM LOGICAL QCOND C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GENCHI,GENNCH EXTERNAL GENCHI,GENNCH C .. C .. EXECUTABLE STATEMENTS .. QCOND = DFN .LE. 1.0D0 .OR. DFD .LE. 0.0D0 .OR. XNONC .LT. 0.0D0 IF (.NOT. (QCOND)) GO TO 10 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 STOP 'DEGREES OF FREEDOM OR NONCENT PARAM OUR OF RANGE IN GENNF' 10 XNUM = GENNCH(DFN,XNONC)/DFN C GENNF = ( GENNCH( DFN, XNONC ) / DFN ) / ( GENCHI( DFD ) / DFD ) XDEN = GENCHI(DFD)/DFD IF (.NOT. (XDEN.LE. (1.0D-38*XNUM))) GO TO 20 WRITE (*,*) ' GENNF - GENERATED NUMBERS WOULD CAUSE OVERFLOW' WRITE (*,*) ' NUMERATOR ',XNUM,' DENOMINATOR ',XDEN WRITE (*,*) ' GENNF RETURNING 1.0E38' GENNF = 1.0D38 GO TO 30 20 GENNF = XNUM/XDEN 30 RETURN END C------------------------------------------------------------------------ DOUBLE PRECISION FUNCTION GENNOR(AV,SD) 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 .. SCALAR ARGUMENTS .. DOUBLE PRECISION AV,SD C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION SNORM EXTERNAL SNORM C .. C .. EXECUTABLE STATEMENTS .. GENNOR = SD*SNORM() + AV RETURN END C-------------------------------------------------------------------------- SUBROUTINE GENPRM(IARRAY,LARRAY) 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 .. SCALAR ARGUMENTS .. INTEGER LARRAY C .. C .. ARRAY ARGUMENTS .. INTEGER IARRAY(LARRAY) C .. C .. LOCAL SCALARS .. INTEGER I,ITMP,IWHICH C .. C .. EXTERNAL FUNCTIONS .. INTEGER IGNUIN EXTERNAL IGNUIN C .. C .. EXECUTABLE STATEMENTS .. DO 10,I = 1,LARRAY IWHICH = IGNUIN(I,LARRAY) ITMP = IARRAY(IWHICH) IARRAY(IWHICH) = IARRAY(I) IARRAY(I) = ITMP 10 CONTINUE RETURN END C--------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION GENUNF(LOW,HIGH) C********************************************************************** C C DOUBLE PRECISION FUNCTION GENUNF( LOW, HIGH ) C C GENERATE UNIFORM DOUBLE PRECISION BETWEEN LOW AND HIGH C C C FUNCTION C C C GENERATES A DOUBLE PRECISION UNIFORMLY DISTRIBUTED BETWEEN LOW AND HIGH. C C C ARGUMENTS C C C LOW --> LOW BOUND (EXCLUSIVE) ON REAL VALUE TO BE GENERATED C DOUBLE PRECISION LOW C C HIGH --> HIGH BOUND (EXCLUSIVE) ON REAL VALUE TO BE GENERATED C DOUBLE PRECISION HIGH C C********************************************************************** C .. SCALAR ARGUMENTS .. DOUBLE PRECISION HIGH,LOW C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION RANF EXTERNAL RANF C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (LOW.GT.HIGH)) GO TO 10 WRITE (*,*) 'LOW > HIGH IN GENUNF: LOW ',LOW,' HIGH: ',HIGH WRITE (*,*) 'ABORT' STOP 'LOW > HIGH IN GENUNF - ABORT' 10 GENUNF = LOW + (HIGH-LOW)*RANF() RETURN END C----------------------------------------------------------------------------- SUBROUTINE GETCGN(G) INTEGER G 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 INTEGER CURNTG,NUMG SAVE CURNTG PARAMETER (NUMG=32) DATA CURNTG/1/ C G = CURNTG RETURN ENTRY SETCGN(G) C********************************************************************** C C SUBROUTINE SETCGN( G ) C SET GENERATOR C C SETS THE CURRENT GENERATOR TO G. ALL REFERENCES TO A GENERAT 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 ABORT IF GENERATOR NUMBER OUT OF RANGE C 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!' STOP ' GENERATOR NUMBER OUT OF RANGE IN SETCGN' 10 CURNTG = G RETURN END C------------------------------------------------------------------------ SUBROUTINE GETSD(ISEED1,ISEED2) C********************************************************************** C C SUBROUTINE GETSD(G,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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. INTEGER ISEED1,ISEED2 C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNIN EXTERNAL QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C ABORT UNLESS RANDOM NUMBER GENERATOR INITIALIZED IF (QRGNIN()) GO TO 10 WRITE (*,*) ' GETSD CALLED BEFORE RANDOM NUMBER GENERATOR ', + ' INITIALIZED -- ABORT!' STOP ' GETSD CALLED BEFORE RANDOM NUMBER GENERATOR INITIALIZED' 10 CALL GETCGN(G) ISEED1 = CG1(G) ISEED2 = CG2(G) RETURN END C------------------------------------------------------------------------- INTEGER FUNCTION IGNBIN(N,PP) 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 DOUBLE PRECISION 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 SUBROUTINE BTPEC(N,PP,ISEED,JX) C C BINOMIAL RANDOM VARIATE GENERATOR C MEAN .LT. 30 -- INVERSE CDF C MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA C FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE C (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE C THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS. C C BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL. C BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE C RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE C USABLE ALGORITHM. C REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER, C "BINOMIAL RANDOM VARIATE GENERATION," C COMMUNICATIONS OF THE ACM, FORTHCOMING C WRITTEN: SEPTEMBER 1980. C LAST REVISED: MAY 1985, JULY 1987 C REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER C GENERATOR C ARGUMENTS C C N : NUMBER OF BERNOULLI TRIALS (INPUT) C PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT) C ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT) C JX: RANDOMLY GENERATED OBSERVATION (OUTPUT) C C VARIABLES C PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC C NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC C XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC C C P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC C FFM: TEMPORARY VARIABLE EQUAL TO XNP + P C M: INTEGER VALUE OF THE CURRENT MODE C FM: FLOATING POINT VALUE OF THE CURRENT MODE C XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS C P1: AREA OF THE TRIANGLE C C: HEIGHT OF THE PARALLELOGRAMS C XM: CENTER OF THE TRIANGLE C XL: LEFT END OF THE TRIANGLE C XR: RIGHT END OF THE TRIANGLE C AL: TEMPORARY VARIABLE C XLL: RATE FOR THE LEFT EXPONENTIAL TAIL C XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL C P2: AREA OF THE PARALLELOGRAMS C P3: AREA OF THE LEFT EXPONENTIAL TAIL C P4: AREA OF THE RIGHT EXPONENTIAL TAIL C U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE C FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE C FROM THE REGION C V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE C (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR C REJECT THE CANDIDATE VALUE C IX: INTEGER CANDIDATE VALUE C X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC C AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC C K: ABSOLUTE VALUE OF (IX-M) C F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE C ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL C ALSO USED IN THE INVERSE TRANSFORMATION C R: THE RATIO P/Q C G: CONSTANT USED IN CALCULATION OF PROBABILITY C MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION C OF F WHEN IX IS GREATER THAN M C IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT C CALCULATION OF F WHEN IX IS LESS THAN M C I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE C AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND C YNORM: LOGARITHM OF NORMAL BOUND C ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V C C X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE C USED IN THE FINAL ACCEPT/REJECT TEST C C QN: PROBABILITY OF NO SUCCESS IN N TRIALS C C REMARK C IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD C SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME C COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS C ARE NOT INVOLVED. C C ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE C GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE C TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR C C********************************************************************** C C C C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY C C .. C .. SCALAR ARGUMENTS .. DOUBLE PRECISION PP INTEGER N C .. C .. LOCAL SCALARS .. DOUBLE PRECISION AL,ALV,AMAXP,C,F,F1,F2,FFM,FM,G,P,P1,P2,P3, +P4,PSAVE,Q,QN,R,U,V,W,W2,X,X1,X2,XL,XLL,XLR,XM,XNP,XNPQ,XR, +YNORM,Z,Z2 INTEGER I,IX,IX1,K,M,MP,NSAVE C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION RANF EXTERNAL RANF C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DLOG,DMIN1,IABS,IDINT,DSQRT C .. C .. DATA STATEMENTS .. DATA PSAVE,NSAVE/-1.,-1/ C .. C .. EXECUTABLE STATEMENTS .. IF (PP.NE.PSAVE) GO TO 10 IF (N.NE.NSAVE) GO TO 20 IF (XNP-30.) 150,30,30 C C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE C 10 PSAVE = PP P = DMIN1(PSAVE,1.-PSAVE) Q = 1.0D0 - P 20 XNP = N*P NSAVE = N IF (XNP.LT.30.) GO TO 140 FFM = XNP + P M = FFM FM = M XNPQ = XNP*Q P1 = IDINT(2.195D0*DSQRT(XNPQ)-4.6D0*Q) + 0.5D0 XM = FM + 0.5D0 XL = XM - P1 XR = XM + P1 C = 0.134D0 + 20.5D0/ (15.3D0+FM) AL = (FFM-XL)/ (FFM-XL*P) XLL = AL* (1.0D0+.50D0*AL) AL = (XR-FFM)/ (XR*Q) XLR = AL* (1.0D0+.5D0*AL) P2 = P1* (1.0D0+C+C) P3 = P2 + C/XLL P4 = P3 + C/XLR C WRITE(6,100) N,P,P1,P2,P3,P4,XL,XR,XM,FM C 100 FORMAT(I15,4F18.7/5F18.7) C C*****GENERATE VARIATE C 30 U = RANF()*P4 V = RANF() C C TRIANGULAR REGION C IF (U.GT.P1) GO TO 40 IX = XM - P1*V + U GO TO 170 C C PARALLELOGRAM REGION C 40 IF (U.GT.P2) GO TO 50 X = XL + (U-P1)/C V = V*C + 1. - DABS(XM-X)/P1 IF (V.GT.1.0D0 .OR. V.LE.0.0D0) GO TO 30 IX = X GO TO 70 C C LEFT TAIL C 50 IF (U.GT.P3) GO TO 60 IX = XL + DLOG(V)/XLL IF (IX.LT.0) GO TO 30 V = V* (U-P2)*XLL GO TO 70 C C RIGHT TAIL C 60 IX = XR - DLOG(V)/XLR IF (IX.GT.N) GO TO 30 V = V* (U-P3)*XLR C C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST C 70 K = IABS(IX-M) IF (K.GT.20 .AND. K.LT.XNPQ/2-1) GO TO 130 C C EXPLICIT EVALUATION C F = 1.0 R = P/Q G = (N+1)*R IF (M-IX) 80,120,100 80 MP = M + 1 DO 90 I = MP,IX F = F* (G/I-R) 90 CONTINUE GO TO 120 100 IX1 = IX + 1 DO 110 I = IX1,M F = F/ (G/I-R) 110 CONTINUE 120 IF (V-F) 170,170,30 C C SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) C 130 AMAXP = (K/XNPQ)* ((K* (K/3.0D0+.6250D0)+ +.1666666666666D0)/XNPQ+.5D0) YNORM = -K*K/ (2.0D0*XNPQ) ALV = DLOG(V) IF (ALV.LT.YNORM-AMAXP) GO TO 170 IF (ALV.GT.YNORM+AMAXP) GO TO 30 C C STIRLING'S FORMULA TO MACHINE ACCURACY FOR C THE FINAL ACCEPTANCE/REJECTION TEST C X1 = IX + 1 F1 = FM + 1.0D0 Z = N + 1 - FM W = N - IX + 1.0D0 Z2 = Z*Z X2 = X1*X1 F2 = F1*F1 W2 = W*W IF (ALV- (XM*DLOG(F1/X1)+ (N-M+.5)*DLOG(Z/W)+ (IX- +M)*DLOG(W*P/ (X1*Q))+ (13860.0D0- (462.0D0- (132.0D0- (99.0D0- +140.D0/F2)/F2)/F2)/F2)/F1/166320.0D0+ (13860.0D0- (462.0D0- +(132.0D0- (99.0D0-140.0D0/Z2)/Z2)/Z2)/Z2)/Z/166320.0D0+ +(13860.0D0- (462.0D0- (132.0D0- (99.0D0-140.0D0/X2)/X2)/X2)/X2) +/X1/166320.0D0+ (13860.0D0- (462.0D0- (132.0D0- (99.0D0- + 140.0D0/W2)/W2)/W2)/W2)/W/166320.0D0)) 170,170,30 C C INVERSE CDF LOGIC FOR MEAN LESS THAN 30 C 140 QN = Q**N R = P/Q G = R* (N+1) 150 IX = 0 F = QN U = RANF() 160 IF (U.LT.F) GO TO 170 IF (IX.GT.110) GO TO 150 U = U - F IX = IX + 1 F = F* (G/IX-R) GO TO 160 170 IF (PSAVE.GT.0.5D0) IX = N - IX IGNBIN = IX RETURN END C------------------------------------------------------------------- INTEGER FUNCTION IGNLGI() 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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER CURNTG,K,S1,S2,Z LOGICAL QQSSD C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNIN EXTERNAL QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN,INRGCM,RGNQSD,SETALL C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C C IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO. C IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO C THIS ROUTINE 2) A CALL TO SETALL. C IF (.NOT. (QRGNIN())) CALL INRGCM() CALL RGNQSD(QQSSD) IF (.NOT. (QQSSD)) CALL SETALL(1234567890,123456789) C C GET CURRENT GENERATOR C CALL GETCGN(CURNTG) S1 = CG1(CURNTG) S2 = CG2(CURNTG) K = S1/53668 S1 = A1* (S1-K*53668) - K*12211 IF (S1.LT.0) S1 = S1 + M1 K = S2/52774 S2 = A2* (S2-K*52774) - K*3791 IF (S2.LT.0) S2 = S2 + M2 CG1(CURNTG) = S1 CG2(CURNTG) = S2 Z = S1 - S2 IF (Z.LT.1) Z = Z + M1 - 1 IF (QANTI(CURNTG)) Z = M1 - Z IGNLGI = Z RETURN END C---------------------------------------------------------------------------- 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 C P --> THE PROBABILITY OF AN EVENT DURING A BERNOULLI TRIAL. C DOUBLE PRECISION P 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 .. DOUBLE PRECISION P INTEGER N C .. C .. LOCAL SCALARS .. DOUBLE PRECISION Y,A,R C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION GENGAM INTEGER IGNPOI EXTERNAL GENGAM,IGNPOI C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DBLE C .. C .. EXECUTABLE STATEMENTS .. C CHECK ARGUMENTS IF (N.LT.0) STOP 'N < 0 IN IGNNBN' IF (P.LE.0.0D0) STOP 'P <= 0 IN IGNNBN' IF (P.GE.1.0D0) STOP 'P >= 1 IN IGNNBN' C GENERATE Y, A RANDOM GAMMA (N,(1-P)/P) VARIABLE R = DBLE(N) A = P/ (1.0D0-P) Y = GENGAM(A,R) C GENERATE A RANDOM POISSON(Y) VARIABLE IGNNBN = IGNPOI(Y) RETURN END C----------------------------------------------------------------------- INTEGER FUNCTION IGNPOI(MU) 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 DOUBLE PRECISION AV C C GENEXP <-- THE RANDOM DEVIATE. C DOUBLE PRECISION 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**********************************************************************C C C C C C P O I S S O N DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C COMPUTER GENERATION OF POISSON DEVIATES C C FROM MODIFIED NORMAL DISTRIBUTIONS. C C ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C C C C (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE) C C C C**********************************************************************C C C INTEGER FUNCTION IGNPOI(IR,MU) C C INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR C MU=MEAN MU OF THE POISSON DISTRIBUTION C OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION C C C C MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR 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 C C C SEPARATION OF CASES A AND B C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION MU C .. C .. LOCAL SCALARS .. DOUBLE PRECISION 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 .. C .. LOCAL ARRAYS .. DOUBLE PRECISION FACT(10),PP(35) C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION RANF,SEXPO,SNORM EXTERNAL RANF,SEXPO,SNORM C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DLOG,DEXP,FLOAT,IFIX,MAX0,MIN0,DSIGN,DSQRT C .. C .. DATA STATEMENTS .. DATA MUPREV,MUOLD/0.,0./ DATA A0,A1,A2,A3,A4,A5,A6,A7/-.5D0,.3333333D0, +-.2500068D0,.2000118D0,-.1661269,.1421878,-.1384794, +.1250060/ DATA FACT/1.0D0,1.0D0,2.0D0,6.0D0,24.0D0,120.0D0,720.0D0, + 5040.0D0,40320.0D0,362880.0D0/ 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 MUPREV = MU S = DSQRT(MU) D = 6.0D0*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 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . C L = IFIX(MU-1.1484D0) C C STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE C 10 G = MU + S*SNORM() IF (G.LT.0.0D0) GO TO 20 IGNPOI = IFIX(G) C C STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH C IF (IGNPOI.GE.L) RETURN C C STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U C FK = FLOAT(IGNPOI) DIFMUK = MU - FK U = RANF() IF (D*U.GE.DIFMUK*DIFMUK*DIFMUK) RETURN C C STEP P. PREPARATIONS FOR STEPS Q AND H. C (RECALCULATIONS OF PARAMETERS IF NECESSARY) C .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. C THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE C APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. C C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. C 20 IF (MU.EQ.MUOLD) GO TO 30 MUOLD = MU OMEGA = .3989423D0/S B1 = .4166667D-1/MU B2 = .3D0*B1*B1 C3 = .1428571D0*B1*B2 C2 = B2 - 15.0D0*C3 C1 = B1 - 6.0D0*B2 + 45.0D0*C3 C0 = 1.0D0 - B1 + 3.0D0*B2 - 15.0D0*C3 C = .1069D0/MU 30 IF (G.LT.0.0D0) GO TO 50 C C 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) C KFLAG = 0 GO TO 70 C C STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) C 40 IF (FY-U*FY.LE.PY*DEXP(PX-FX)) RETURN C C STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL C DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' C (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) C 50 E = SEXPO() U = RANF() U = U + U - 1.0D0 T = 1.8D0 + DSIGN(E,U) IF (T.LE. (-.6744D0)) GO TO 50 IGNPOI = IFIX(MU+S*T) FK = FLOAT(IGNPOI) DIFMUK = MU - FK C C 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) C KFLAG = 1 GO TO 70 C C STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) C 60 IF (C*DABS(U).GT.PY*DEXP(PX+E)-FY*DEXP(FX+E)) GO TO 50 RETURN C C STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY. C CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT C 70 IF (IGNPOI.GE.10) GO TO 80 PX = -MU PY = MU**IGNPOI/FACT(IGNPOI+1) GO TO 110 C C CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION C A0-A7 FOR ACCURACY WHEN ADVISABLE C .8333333E-1=1./12. .3989423=(2*PI)**(-.5) C 80 DEL = .8333333D-1/FK DEL = DEL - 4.8D0*DEL*DEL*DEL V = DIFMUK/FK IF (DABS(V).LE.0.25D0) GO TO 90 PX = FK*DLOG(1.0D0+V) - DIFMUK - DEL GO TO 100 90 PX = FK*V*V* (((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V+A0) - + DEL 100 PY = .3989423D0/DSQRT(FK) 110 X = (0.5D0-DIFMUK)/S XX = X*X FX = -0.5D0*XX FY = OMEGA* (((C3*XX+C2)*XX+C1)*XX+C0) IF (KFLAG) 40,40,60 C C C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) C 120 MUPREV = 0.0D0 IF (MU.EQ.MUOLD) GO TO 130 MUOLD = MU M = MAX0(1,IFIX(MU)) L = 0 P = DEXP(-MU) Q = P P0 = P C C STEP U. UNIFORM SAMPLE FOR INVERSION METHOD C 130 U = RANF() IGNPOI = 0 IF (U.LE.P0) RETURN C C STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE C PP-TABLE OF CUMULATIVE POISSON PROBABILITIES C (0.458=PP(9) FOR MU=10) C IF (L.EQ.0) GO TO 150 J = 1 IF (U.GT.0.458D0) J = MIN0(L,M) DO 140 K = J,L IF (U.LE.PP(K)) GO TO 180 140 CONTINUE IF (L.EQ.35) GO TO 130 C C STEP C. CREATION OF NEW POISSON PROBABILITIES P C AND THEIR CUMULATIVES Q=PP(K) C 150 L = L + 1 DO 160 K = L,35 P = P*MU/FLOAT(K) Q = Q + P PP(K) = Q IF (U.LE.Q) GO TO 170 160 CONTINUE L = 35 GO TO 130 170 L = K 180 IGNPOI = K RETURN END C------------------------------------------------------------------- INTEGER FUNCTION IGNUIN(LOW,HIGH) 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 IGNLGI GENERATES INTEGERS BETWEEN 1 AND 2147483562 C MAXNUM IS 1 LESS THAN MAXIMUM GENERABLE VALUE C .. PARAMETERS .. INTEGER MAXNUM PARAMETER (MAXNUM=2147483561) CHARACTER*(*) ERR1,ERR2 PARAMETER (ERR1='LOW > HIGH IN IGNUIN', + ERR2=' ( HIGH - LOW ) > 2,147,483,561 IN IGNUIN') C .. C .. SCALAR ARGUMENTS .. INTEGER HIGH,LOW C .. C .. LOCAL SCALARS .. INTEGER ERR,IGN,MAXNOW,RANGE,RANP1 C .. C .. EXTERNAL FUNCTIONS .. INTEGER IGNLGI EXTERNAL IGNLGI C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC MOD C .. C .. EXECUTABLE STATEMENTS .. IF (.NOT. (LOW.GT.HIGH)) GO TO 10 ERR = 1 C ABORT-PROGRAM GO TO 80 10 RANGE = HIGH - LOW IF (.NOT. (RANGE.GT.MAXNUM)) GO TO 20 ERR = 2 C ABORT-PROGRAM GO TO 80 20 IF (.NOT. (LOW.EQ.HIGH)) GO TO 30 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 30 RANP1 = RANGE + 1 MAXNOW = (MAXNUM/RANP1)*RANP1 40 IGN = IGNLGI() - 1 IF (.NOT. (IGN.LE.MAXNOW)) GO TO 50 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 C TO ABORT-PROGRAM 90 WRITE (*,*) ERR2 100 WRITE (*,*) ' LOW: ',LOW,' HIGH: ',HIGH WRITE (*,*) ' ABORT ON FATAL ERROR' IF (.NOT. (ERR.EQ.1)) GO TO 110 STOP 'LOW > HIGH IN IGNUIN' GO TO 120 110 STOP ' ( HIGH - LOW ) > 2,147,483,561 IN IGNUIN' 120 END C------------------------------------------------------------------- SUBROUTINE INITGN(ISDTYP) C********************************************************************** C C SUBROUTINE INITGN(ISDTYP) C INIT-IALIZE CURRENT G-E-N-ERATOR C C REINITIALIZES THE STATE OF THE CURRENT GENERATOR C C THIS IS A TRANSCRIPTION FROM PASCAL TO FORTRAN OF ROUTINE C INIT_GENERATOR FROM THE PAPER C C L'ECUYER, P. AND COTE, S. "IMPLEMENTING A RANDOM NUMBER PACKAGE C WITH SPLITTING FACILITIES." ACM TRANSACTIONS ON MATHEMATICAL C SOFTWARE, 17:98-111 (1991) C C C ARGUMENTS C C C ISDTYP -> THE STATE TO WHICH THE GENERATOR IS TO BE SET C C ISDTYP = -1 => SETS THE SEEDS TO THEIR INITIAL VALUE C ISDTYP = 0 => SETS THE SEEDS TO THE FIRST VALUE OF C THE CURRENT BLOCK C ISDTYP = 1 => SETS THE SEEDS TO THE FIRST VALUE OF C THE NEXT BLOCK C C INTEGER ISDTYP C C********************************************************************** C .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. INTEGER ISDTYP C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNIN INTEGER MLTMOD EXTERNAL QRGNIN,MLTMOD C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C ABORT UNLESS RANDOM NUMBER GENERATOR INITIALIZED IF (QRGNIN()) GO TO 10 WRITE (*,*) ' INITGN CALLED BEFORE RANDOM NUMBER GENERATOR ', + ' INITIALIZED -- ABORT!' STOP ' INITGN CALLED BEFORE RANDOM NUMBER GENERATOR INITIALIZED' 10 CALL GETCGN(G) IF ((-1).NE. (ISDTYP)) GO TO 20 LG1(G) = IG1(G) LG2(G) = IG2(G) GO TO 50 20 IF ((0).NE. (ISDTYP)) GO TO 30 CONTINUE GO TO 50 C DO NOTHING 30 IF ((1).NE. (ISDTYP)) GO TO 40 LG1(G) = MLTMOD(A1W,LG1(G),M1) LG2(G) = MLTMOD(A2W,LG2(G),M2) GO TO 50 40 STOP 'ISDTYP NOT IN RANGE' 50 CG1(G) = LG1(G) CG2(G) = LG2(G) RETURN END C----------------------------------------------------------------------------- SUBROUTINE INRGCM() 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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER I LOGICAL QDUM C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNSN EXTERNAL QRGNSN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C V=20; W=30; C C A1W = MOD(A1**(2**W),M1) A2W = MOD(A2**(2**W),M2) C A1VW = MOD(A1**(2**(V+W)),M1) A2VW = MOD(A2**(2**(V+W)),M2) C C IF V OR W IS CHANGED A1W, A2W, A1VW, AND A2VW NEED TO BE RECOMPUTED. C AN EFFICIENT WAY TO PRECOMPUTE A**(2*J) MOD M IS TO START WITH C A AND SQUARE IT J TIMES MODULO M USING THE FUNCTION MLTMOD. C M1 = 2147483563 M2 = 2147483399 A1 = 40014 A2 = 40692 A1W = 1033780774 A2W = 1494757890 A1VW = 2082007225 A2VW = 784306273 DO 10,I = 1,NUMG QANTI(I) = .FALSE. 10 CONTINUE C C TELL THE WORLD THAT COMMON HAS BEEN INITIALIZED C QDUM = QRGNSN(.TRUE.) RETURN END C----------------------------------------------------------------------- INTEGER FUNCTION LENNOB(STRING) IMPLICIT INTEGER (A-P,R-Z),LOGICAL (Q) C********************************************************************** C C INTEGER FUNCTION LENNOB( STRING ) C LENGTH NOT COUNTING TRAILING BLANKS C C C FUNCTION C C C RETURNS THE LENGTH OF STRING UP TO AND INCLUDING THE LAST C NON-BLANK CHARACTER. C C C ARGUMENTS C C C STRING --> STRING WHOSE LENGTH NOT COUNTING TRAILING BLANKS C IS RETURNED. C C********************************************************************** CHARACTER*(*) STRING END = LEN(STRING) DO 20,I = END,1,-1 IF (.NOT. (STRING(I:I).NE.' ')) GO TO 10 LENNOB = I RETURN 10 CONTINUE 20 CONTINUE LENNOB = 0 RETURN END C-------------------------------------------------------------------------- INTEGER FUNCTION MLTMOD(A,S,M) 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 .. PARAMETERS .. INTEGER H PARAMETER (H=32768) C .. C .. SCALAR ARGUMENTS .. INTEGER A,M,S C .. C .. LOCAL SCALARS .. INTEGER A0,A1,K,P,Q,QH,RH C .. C .. EXECUTABLE STATEMENTS .. C C H = 2**((B-2)/2) WHERE B = 32 BECAUSE WE ARE USING A 32 BIT C MACHINE. ON A DIFFERENT MACHINE RECOMPUTE H C IF (.NOT. (A.LE.0.OR.A.GE.M.OR.S.LE.0.OR.S.GE.M)) GO TO 10 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' STOP ' A, M, S OUT OF ORDER IN MLTMOD - ABORT!' 10 IF (.NOT. (A.LT.H)) GO TO 20 A0 = A P = 0 GO TO 120 20 A1 = A/H A0 = A - H*A1 QH = M/H RH = M - H*QH IF (.NOT. (A1.GE.H)) GO TO 50 A1 = A1 - H K = S/QH P = H* (S-K*QH) - K*RH 30 IF (.NOT. (P.LT.0)) GO TO 40 P = P + M GO TO 30 40 GO TO 60 50 P = 0 C C P = (A2*S*H)MOD M C 60 IF (.NOT. (A1.NE.0)) GO TO 90 Q = M/A1 K = S/Q P = P - K* (M-A1*Q) IF (P.GT.0) P = P - M P = P + A1* (S-K*Q) 70 IF (.NOT. (P.LT.0)) GO TO 80 P = P + M GO TO 70 80 CONTINUE 90 K = P/QH C C P = ((A2*H + A1)*S)MOD M C P = H* (P-K*QH) - K*RH 100 IF (.NOT. (P.LT.0)) GO TO 110 P = P + M GO TO 100 110 CONTINUE 120 IF (.NOT. (A0.NE.0)) GO TO 150 C C P = ((A2*H + A1)*H*S)MOD M C Q = M/A0 K = S/Q P = P - K* (M-A0*Q) IF (P.GT.0) P = P - M P = P + A0* (S-K*Q) 130 IF (.NOT. (P.LT.0)) GO TO 140 P = P + M GO TO 130 140 CONTINUE 150 MLTMOD = P C RETURN END C-------------------------------------------------------------------- SUBROUTINE PHRTSD(PHRASE,SEED1,SEED2) 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 .. PARAMETERS .. CHARACTER*(*) TABLE PARAMETER (TABLE='ABCDEFGHIJKLMNOPQRSTUVWXYZ'// + 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'//'0123456789'// + '!@#$%^&*()_+[];:''"<>?,./') INTEGER TWOP30 PARAMETER (TWOP30=1073741824) C .. C .. SCALAR ARGUMENTS .. INTEGER SEED1,SEED2 CHARACTER PHRASE* (*) C .. C .. LOCAL SCALARS .. INTEGER I,ICHR,J,LPHR C .. C .. LOCAL ARRAYS .. INTEGER SHIFT(0:4),VALUES(5) C .. C .. EXTERNAL FUNCTIONS .. INTEGER LENNOB EXTERNAL LENNOB C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC INDEX,MOD C .. C .. DATA STATEMENTS .. DATA SHIFT/1,64,4096,262144,16777216/ C .. C .. EXECUTABLE STATEMENTS .. SEED1 = 1234567890 SEED2 = 123456789 LPHR = LENNOB(PHRASE) IF (LPHR.LT.1) RETURN DO 30,I = 1,LPHR ICHR = MOD(INDEX(TABLE,PHRASE(I:I)),64) IF (ICHR.EQ.0) ICHR = 63 DO 10,J = 1,5 VALUES(J) = ICHR - J IF (VALUES(J).LT.1) VALUES(J) = VALUES(J) + 63 10 CONTINUE DO 20,J = 1,5 SEED1 = MOD(SEED1+SHIFT(J-1)*VALUES(J),TWOP30) SEED2 = MOD(SEED2+SHIFT(J-1)*VALUES(6-J),TWOP30) 20 CONTINUE 30 CONTINUE RETURN END C------------------------------------------------------------------------- LOGICAL FUNCTION QRGNIN() C********************************************************************** C C LOGICAL FUNCTION QRGNIN() C Q RANDOM GENERATORS INITIALIZED? C C A TRIVIAL ROUTINE TO DETERMINE WHETHER OR NOT THE RANDOM C NUMBER GENERATOR HAS BEEN INITIALIZED. RETURNS .TRUE. IF C IT HAS, ELSE .FALSE. C C********************************************************************** C .. SCALAR ARGUMENTS .. LOGICAL QVALUE C .. C .. LOCAL SCALARS .. LOGICAL QINIT C .. C .. ENTRY POINTS .. LOGICAL QRGNSN C .. C .. SAVE STATEMENT .. SAVE QINIT C .. C .. DATA STATEMENTS .. DATA QINIT/.FALSE./ C .. C .. EXECUTABLE STATEMENTS .. QRGNIN = QINIT RETURN ENTRY QRGNSN(QVALUE) C********************************************************************** C C LOGICAL FUNCTION QRGNSN( QVALUE ) C Q RANDOM GENERATORS SET WHETHER INITIALIZED C C SETS STATE OF WHETHER RANDOM NUMBER GENERATOR IS INITIALIZED C TO QVALUE. C C THIS ROUTINE IS ACTUALLY AN ENTRY IN QRGNIN, HENCE IT IS A C LOGICAL FUNCTION. IT RETURNS THE (MEANINGLESS) VALUE .TRUE. C C********************************************************************** QINIT = QVALUE QRGNSN = .TRUE. RETURN END C---------------------------------------------------------------------- DOUBLE PRECISION FUNCTION RANF() C********************************************************************** C C DOUBLE PRECISION 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 .. EXTERNAL FUNCTIONS .. INTEGER IGNLGI EXTERNAL IGNLGI C .. C .. EXECUTABLE STATEMENTS .. C C 4.656613057E-10 IS 1/M1 M1 IS SET IN A DATA STATEMENT IN IGNLGI C AND IS CURRENTLY 2147483563. IF M1 CHANGES, CHANGE THIS ALSO. C RANF = IGNLGI()*4.656613057D-10 RETURN END C-------------------------------------------------------------------------- SUBROUTINE SETALL(ISEED1,ISEED2) 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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. INTEGER ISEED1,ISEED2 LOGICAL QSSD C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G,OCGN LOGICAL QQSSD C .. C .. EXTERNAL FUNCTIONS .. INTEGER MLTMOD LOGICAL QRGNIN EXTERNAL MLTMOD,QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN,INITGN,INRGCM,SETCGN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/,QQSSD C .. C .. DATA STATEMENTS .. DATA QQSSD/.FALSE./ C .. C .. EXECUTABLE STATEMENTS .. C C TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE C HAS BEEN CALLED. C QQSSD = .TRUE. CALL GETCGN(OCGN) C C INITIALIZE COMMON BLOCK IF NECESSARY C IF (.NOT. (QRGNIN())) CALL INRGCM() IG1(1) = ISEED1 IG2(1) = ISEED2 CALL INITGN(-1) DO 10,G = 2,NUMG IG1(G) = MLTMOD(A1VW,IG1(G-1),M1) IG2(G) = MLTMOD(A2VW,IG2(G-1),M2) CALL SETCGN(G) CALL INITGN(-1) 10 CONTINUE CALL SETCGN(OCGN) RETURN ENTRY RGNQSD(QSSD) C********************************************************************** C C SUBROUTINE RGNQSD C RANDOM NUMBER GENERATOR QUERY SEED SET? C C RETURNS (LOGICAL) QSSD AS .TRUE. IF SETALL HAS BEEN INVOKED, C OTHERWISE RETURNS .FALSE. C C********************************************************************** QSSD = QQSSD RETURN END C------------------------------------------------------------------------- SUBROUTINE SETANT(QVALUE) 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 .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. LOGICAL QVALUE C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNIN EXTERNAL QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C ABORT UNLESS RANDOM NUMBER GENERATOR INITIALIZED IF (QRGNIN()) GO TO 10 WRITE (*,*) ' SETANT CALLED BEFORE RANDOM NUMBER GENERATOR ', + ' INITIALIZED -- ABORT!' STOP ' SETANT CALLED BEFORE RANDOM NUMBER GENERATOR INITIALIZED' 10 CALL GETCGN(G) QANTI(G) = QVALUE RETURN END C------------------------------------------------------------------------- SUBROUTINE SETGMN(MEANV,CHOL,P,PARM) 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 DOUBLE PRECISION PARM(P*(P+3)/2 + 1) C C********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER P C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION CHOL(P,P),MEANV(P),PARM(P* (P+3)/2+1) C .. C .. LOCAL SCALARS .. INTEGER I,ICOUNT,INFO,J C .. C .. EXTERNAL SUBROUTINES .. C .. C .. EXECUTABLE STATEMENTS .. C C C TEST THE INPUT C IF (.NOT. (P.LE.0)) GO TO 10 WRITE (*,*) 'P NONPOSITIVE IN SETGMN' WRITE (*,*) 'VALUE OF P: ',P STOP 'P NONPOSITIVE IN SETGMN' 10 PARM(1) = P C C PUT P AND MEANV INTO PARM C DO 20,I = 2,P + 1 PARM(I) = MEANV(I-1) 20 CONTINUE 30 ICOUNT = P + 1 C C PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM C COVM(1,1) = PARM(P+2) C COVM(1,2) = PARM(P+3) C : C COVM(1,P) = PARM(2P+1) C COVM(2,2) = PARM(2P+2) ... C DO 50,I = 1,P DO 40,J = I,P ICOUNT = ICOUNT + 1 PARM(ICOUNT) = CHOL(I,J) 40 CONTINUE 50 CONTINUE RETURN C END C----------------------------------------------------------------------------- SUBROUTINE SETSD(ISEED1,ISEED2) 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 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********************************************************************** C .. PARAMETERS .. INTEGER NUMG PARAMETER (NUMG=32) C .. C .. SCALAR ARGUMENTS .. INTEGER ISEED1,ISEED2 C .. C .. SCALARS IN COMMON .. INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2 C .. C .. ARRAYS IN COMMON .. INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG), + LG2(NUMG) LOGICAL QANTI(NUMG) C .. C .. LOCAL SCALARS .. INTEGER G C .. C .. EXTERNAL FUNCTIONS .. LOGICAL QRGNIN EXTERNAL QRGNIN C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GETCGN,INITGN C .. C .. COMMON BLOCKS .. COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1, + CG2,QANTI C .. C .. SAVE STATEMENT .. SAVE /GLOBE/ C .. C .. EXECUTABLE STATEMENTS .. C ABORT UNLESS RANDOM NUMBER GENERATOR INITIALIZED IF (QRGNIN()) GO TO 10 WRITE (*,*) ' SETSD CALLED BEFORE RANDOM NUMBER GENERATOR ', + ' INITIALIZED -- ABORT!' STOP ' SETSD CALLED BEFORE RANDOM NUMBER GENERATOR INITIALIZED' 10 CALL GETCGN(G) IG1(G) = ISEED1 IG2(G) = ISEED2 CALL INITGN(-1) RETURN END C--------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SEXPO() C**********************************************************************C C C C C C (STANDARD-) E X P O N E N T I A L DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C COMPUTER METHODS FOR SAMPLING FROM THE C C EXPONENTIAL AND NORMAL DISTRIBUTIONS. C C COMM. ACM, 15,10 (OCT. 1972), 873 - 882. C C C C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM C C 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C C C C MODIFIED BY BARRY W. BROWN, FEB 3, 1988 TO USE RANF INSTEAD OF C C SUNIF. THE ARGUMENT IR THUS GOES AWAY. C 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 DATA Q/.6931472D0,.9333737D0,.9888778D0,.9984959D0,.9998293D0, +.9999833D0,.9999986D0,.9999999D0/ C 10 A = 0.0D0 U = RANF() GO TO 30 20 A = A + Q1 30 U = U + U IF (U.LE.1.0D0) GO TO 20 40 U = U - 1.0D0 IF (U.GT.Q1) GO TO 60 50 SEXPO = A + U RETURN 60 I = 1 USTAR = RANF() UMIN = USTAR 70 USTAR = RANF() IF (USTAR.LT.UMIN) UMIN = USTAR 80 I = I + 1 IF (U.GT.Q(I)) GO TO 70 90 SEXPO = A + UMIN*Q1 RETURN END C------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SGAMMA(A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C**********************************************************************C C C C C C (STANDARD-) G A M M A DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C PARAMETER A >= 1.0 ! C C C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C GENERATING GAMMA VARIATES BY A C C MODIFIED REJECTION TECHNIQUE. C C COMM. ACM, 25,1 (JAN. 1982), 47 - 54. C C C C STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER C C (STRAIGHTFORWARD IMPLEMENTATION) C C C C MODIFIED BY BARRY W. BROWN, FEB 3, 1988 TO USE RANF INSTEAD OF C C SUNIF. THE ARGUMENT IR THUS GOES AWAY. C C C C**********************************************************************C C C C PARAMETER 0.0 < A < 1.0 ! C C C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C COMPUTER METHODS FOR SAMPLING FROM GAMMA, C C BETA, POISSON AND BINOMIAL DISTRIBUTIONS. C C COMPUTING, 12 (1974), 223 - 246. C C C C (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER) C C C C**********************************************************************C C C C INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION C OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION C C COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K)) 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 DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/.04166669D0,.02083148D0,.00801191D0, +.00144121D0,-.00007388D0,.00024511D0,.00024240D0/ DATA A1,A2,A3,A4,A5,A6,A7/.3333333D0,-.2500030D0,.2000062D0, +-.1662921D0,.1423657D0,-.1367177D0,.1233795D0/ DATA E1,E2,E3,E4,E5/1.0D0,.4999897D0,.1668290D0,.0407753D0, +.0102930D0/ 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.0D0/,AAA/0.0D0/,SQRT32/5.656854D0/ C C SAVE STATEMENTS C SAVE AA,AAA,S2,S,D,Q0,B,SI,C C IF (A.EQ.AA) GO TO 10 IF (A.LT.1.0) GO TO 140 C C STEP 1: RECALCULATIONS OF S2,S,D IF A HAS CHANGED C AA = A S2 = A - 0.5D0 S = DSQRT(S2) D = SQRT32 - 12.0D0*S C C STEP 2: T=STANDARD NORMAL DEVIATE, C X=(S,1/2)-NORMAL DEVIATE. C IMMEDIATE ACCEPTANCE (I) C 10 T = SNORM() X = S + 0.5D0*T SGAMMA = X*X IF (T.GE.0.0D0) RETURN C C STEP 3: U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S) C U = RANF() IF (D*U.LE.T*T*T) RETURN C C STEP 4: RECALCULATIONS OF Q0,B,SI,C IF NECESSARY C IF (A.EQ.AAA) GO TO 40 AAA = A R = 1.0D0/A Q0 = ((((((Q7*R+Q6)*R+Q5)*R+Q4)*R+Q3)*R+Q2)*R+Q1)*R C C APPROXIMATION DEPENDING ON SIZE OF PARAMETER A C THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND C C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS C IF (A.LE.3.686D0) GO TO 30 IF (A.LE.13.022D0) GO TO 20 C C CASE 3: A .GT. 13.022 C B = 1.77D0 SI = .75D0 C = .1515D0/S GO TO 40 C C CASE 2: 3.686 .LT. A .LE. 13.022 C 20 B = 1.654D0 + .0076D0*S2 SI = 1.68D0/S + .275D0 C = .062D0/S + .024D0 GO TO 40 C C CASE 1: A .LE. 3.686 C 30 B = .463D0 + S + .178D0*S2 SI = 1.235D0 C = .195D0/S - .079D0 + .16D0*S C C STEP 5: NO QUOTIENT TEST IF X NOT POSITIVE C 40 IF (X.LE.0.0D0) GO TO 70 C C STEP 6: CALCULATION OF V AND QUOTIENT Q C V = T/ (S+S) IF (DABS(V).LE.0.25D0) GO TO 50 Q = Q0 - S*T + 0.25D0*T*T + (S2+S2)*DLOG(1.0D0+V) GO TO 60 50 Q = Q0 + 0.5D0*T*T* ((((((A7*V+A6)*V+A5)*V+A4)*V+A3) +*V+A2)*V+A1)*V C C STEP 7: QUOTIENT ACCEPTANCE (Q) C 60 IF (DLOG(1.0D0-U).LE.Q) RETURN C C STEP 8: E=STANDARD EXPONENTIAL DEVIATE C U= 0,1 -UNIFORM DEVIATE C T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE C 70 E = SEXPO() U = RANF() U = U + U - 1.0D0 T = B + DSIGN(SI*E,U) IF (.NOT. (U.GE.0.0D0)) GO TO 80 T = B + SI*E GO TO 90 80 T = B - SI*E C C STEP 9: REJECTION IF T .LT. TAU(1) = -.71874483771719 C 90 IF (T.LT. (-.7187449D0)) GO TO 70 C C STEP 10: CALCULATION OF V AND QUOTIENT Q C V = T/ (S+S) IF (DABS(V).LE.0.25D0) GO TO 100 Q = Q0 - S*T + 0.25D0*T*T + (S2+S2)*DLOG(1.0D0+V) GO TO 110 100 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 110 IF (Q.LE.0.0D0) GO TO 70 IF (Q.LE.0.5D0) GO TO 120 W = DEXP(Q) - 1.0D0 GO TO 130 120 W = ((((E5*Q+E4)*Q+E3)*Q+E2)*Q+E1)*Q C C IF T IS REJECTED, SAMPLE AGAIN AT STEP 8 C 130 IF (C*DABS(U).GT.W*DEXP(E-0.5D0*T*T)) GO TO 70 X = S + 0.5D0*T SGAMMA = X*X RETURN C C ALTERNATE METHOD FOR PARAMETERS A BELOW 1 (.3678794=EXP(-1.)) C 140 AA = 0.0D0 B = 1.0D0 + .3678794D0*A 150 P = B*RANF() IF (P.GE.1.0) GO TO 160 SGAMMA = DEXP(DLOG(P)/A) IF (SEXPO().LT.SGAMMA) GO TO 150 RETURN 160 SGAMMA = -DLOG((B-P)/A) IF (SEXPO().LT. (1.0D0-A)*DLOG(SGAMMA)) GO TO 150 RETURN END C---------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SNORM() IMPLICIT DOUBLE PRECISION(A-H,O-Z) C**********************************************************************C C C C C C (STANDARD-) N O R M A L DISTRIBUTION C C C C C C**********************************************************************C C**********************************************************************C C C C FOR DETAILS SEE: C C C C AHRENS, J.H. AND DIETER, U. C C EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM C C SAMPLING FROM THE NORMAL DISTRIBUTION. C C MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. C C C C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' C C (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C C C C MODIFIED BY BARRY W. BROWN, FEB 3, 1988 TO USE RANF INSTEAD OF C C SUNIF. THE ARGUMENT IR THUS GOES AWAY. C 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 DATA A/0.0D0,.3917609D-1,.7841241D-1,.1177699D0,.1573107D0, +.1970991D0,.2372021D0,.2776904D0,.3186394D0,.3601299D0, +.4022501D0,.4450965D0,.4887764D0,.5334097D0,.5791322D0,.6260990D0, +.6744898D0,.7245144D0,.7764218D0,.8305109D0,.8871466D0, +.9467818D0,1.009990D0,1.077516D0,1.150349D0,1.229859D0,1.318011D0, + 1.417797D0,1.534121D0,1.675940D0,1.862732D0,2.153875D0/ DATA D/5*0.0D0,.2636843D0,.2425085D0,.2255674D0,.2116342D0, +.1999243D0,.1899108D0,.1812252D0,.1736014D0,.1668419D0, +.1607967D0,.1553497D0,.1504094D0,.1459026D0,.1417700D0, +.1379632D0,.1344418D0,.1311722D0,.1281260D0,.1252791D0, +.1226109D0,.1201036D0,.1177417D0,.1155119D0,.1134023D0, +.1114027D0,.1095039D0/ DATA T/.7673828D-3,.2306870D-2,.3860618D-2,.5438454D-2, + .7050699D-2,.8708396D-2,.1042357D-1,.1220953D-1,.1408125D-1, + .1605579D-1,.1815290D-1,.2039573D-1,.2281177D-1,.2543407D-1, + .2830296D-1,.3146822D-1,.3499233D-1,.3895483D-1,.4345878D-1, + .4864035D-1,.5468334D-1,.6184222D-1,.7047983D-1,.8113195D-1, + .9462444D-1,.1123001D0,.1364980D0,.1716886D0,.2276241D0, +.3304980D0,.5847031D0/ DATA H/.3920617D-1,.3932705D-1,.3950999D-1,.3975703D-1, + .4007093D-1,.4045533D-1,.4091481D-1,.4145507D-1,.4208311D-1, + .4280748D-1,.4363863D-1,.4458932D-1,.4567523D-1,.4691571D-1, + .4833487D-1,.4996298D-1,.5183859D-1,.5401138D-1,.5654656D-1, + .5953130D-1,.6308489D-1,.6737503D-1,.7264544D-1,.7926471D-1, + .8781922D-1,.9930398D-1,.1155599D0,.1404344D0,.1836142D0, +.2790016D0,.7010474D0/ C 10 U = RANF() S = 0.0D0 IF (U.GT.0.5D0) S = 1.0D0 U = U + U - S 20 U = 32.0D0*U I = IDINT(U) IF (I.EQ.32) I = 31 IF (I.EQ.0) GO TO 100 C C START CENTER C 30 USTAR = U - FLOAT(I) AA = A(I) 40 IF (USTAR.LE.T(I)) GO TO 60 W = (USTAR-T(I))*H(I) C C EXIT (BOTH CASES) C 50 Y = AA + W SNORM = Y IF (S.EQ.1.0D0) SNORM = -Y RETURN C C CENTER CONTINUED C 60 U = RANF() W = U* (A(I+1)-AA) TT = (0.5D0*W+AA)*W GO TO 80 70 TT = U USTAR = RANF() 80 IF (USTAR.GT.TT) GO TO 50 90 U = RANF() IF (USTAR.GE.U) GO TO 70 USTAR = RANF() GO TO 40 C C START TAIL C 100 I = 6 AA = A(32) GO TO 120 110 AA = AA + D(I) I = I + 1 120 U = U + U IF (U.LT.1.0D0) GO TO 110 130 U = U - 1.0D0 140 W = U*D(I) TT = (0.5D0*W+AA)*W GO TO 160 150 TT = U 160 USTAR = RANF() IF (USTAR.GT.TT) GO TO 50 170 U = RANF() IF (USTAR.GE.U) GO TO 150 U = RANF() GO TO 140 END C-------------------------------------------------------------------------- *DECK SPOFA SUBROUTINE SPOFA(A,LDA,N,INFO) INTEGER LDA,N,INFO DOUBLE PRECISION A(LDA,1) C C SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX. C C SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) . C C ON ENTRY C C A REAL(LDA, N) C THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE C DIAGONAL AND UPPER TRIANGLE ARE USED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R C WHERE TRANS(R) IS THE TRANSPOSE. C THE STRICT LOWER TRIANGLE IS UNALTERED. C IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE. C C INFO INTEGER C = 0 FOR NORMAL RETURN. C = K SIGNALS AN ERROR CONDITION. THE LEADING MINOR C OF ORDER K IS NOT POSITIVE DEFINITE. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SDOT C FORTRAN SQRT C C INTERNAL VARIABLES C DOUBLE PRECISION SDOT,T DOUBLE PRECISION S INTEGER J,JM1,K C BEGIN BLOCK WITH ...EXITS TO 40 C C DO 30 J = 1,N INFO = J S = 0.0D0 JM1 = J - 1 IF (JM1.LT.1) GO TO 20 DO 10 K = 1,JM1 T = A(K,J) - SDOT(K-1,A(1,K),1,A(1,J),1) T = T/A(K,K) A(K,J) = T S = S + T*T 10 CONTINUE 20 CONTINUE S = A(J,J) - S C ......EXIT IF (S.LE.0.0D0) GO TO 40 A(J,J) = DSQRT(S) 30 CONTINUE INFO = 0 40 CONTINUE RETURN END C---------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SDOT(N,SX,INCX,SY,INCY) DOUBLE PRECISION SX(1),SY(1),STEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N STEMP = 0.0D0 SDOT = 0.0D0 IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE SDOT = STEMP RETURN 20 M = MOD(N,5) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M STEMP = STEMP + SX(I)*SY(I) 30 CONTINUE IF (N.LT.5) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 STEMP = STEMP + SX(I)*SY(I) + SX(I+1)*SY(I+1) + + SX(I+2)*SY(I+2) + SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4) 50 CONTINUE 60 SDOT = STEMP RETURN END