C ALGORITHM 564 C C A TEST PROBLEM GENERATOR FOR DISCRETE LINEAR L1 APPROXIMATION C PROBLEMS C C BY K.L. HOFFMAN AND D.R. SHIER C C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 6 (DECEMBER 1980), C 615-617 C C MAN 10 C THIS IS A DRIVER ROUTINE TO GENERATE A DATA SET WHICH CAN MAN 20 C BE USED FOR TESTING L1-APPROXIMATION (LEAST ABSOLUTE-DEVIATION MAN 30 C CURVE FITTING) COMPUTER CODES. MAN 40 C MAN 50 C MAN 60 C INPUT ARGUMENTS TO BE READ FROM UNIT 5' MAN 70 C MAN 80 C NOBS -- NUMBER OF OBSERVATIONS MAN 90 C NPAR -- NUMBER OF PARAMETERS MAN 100 C NXGEN -- INDICATES WHETHER SOLUTION X IS SPECIFIED MAN 110 C AS INPUT OR GENERATED MAN 120 C = 0 IF SPECIFIED AS INPUT MAN 130 C = 1 IF RANDOMLY GENERATED FROM NORMAL MAN 140 C = 2 IF RANDOMLY GENERATED FROM UNIFORM MAN 150 C XOPT -- REAL ARRAY OF DIMENSION NPAR - RLOSS WHICH MAN 160 C SPECIFIES THE OPTIMAL SOLUTION IF NXGEN IS MAN 170 C SET EQUAL TO ZERO MAN 180 C XMEAN -- MEAN OF NORMAL DISTRIBUTION IF NXGEN=1, MAN 190 C LOWER LIMIT OF UNIFORM DISTRIBUTION IF MAN 200 C NXGEN=2 MAN 210 C XVAR -- STANDARD DEVIATION OF NORMAL DISTRIBUTION MAN 220 C IF NXGEN=1, UPPER LIMIT OF UNIFORM DISTRI- MAN 230 C BUTION IF NXGEN=2 MAN 240 C NDIST -- ARRAY OF DIMENSION NPAR - RLOSS WHICH MAN 250 C SPECIFIES THE STATISTICAL DISTRIBUTION MAN 260 C FOR EACH COLUMN OF THE DESIGN MATRIX MAN 270 C = 0 IF NORMAL MAN 280 C = 1 IF UNIFORM MAN 290 C = 2 IF ZERO-ONE MAN 300 C CLMEAN -- ARRAY OF DIMENSION NPAR - RLOSS WHICH MAN 310 C SPECIFIES FOR EACH COLUMN THE MEAN (IF MAN 320 C NDIST=0), LOWER LIMIT (IF NDIST=1), OR MAN 330 C PROPORTION P OF ZERO ENTRIES (IF NDIST=2) MAN 340 C CLMVAR -- ARRAY OF DIMENSION NPAR - RLOSS WHICH MAN 350 C SPECIFIES FOR EACH COLUMN THE STANDARD MAN 360 C DEVIATION (IF NDIST=0), THE UPPER LIMIT MAN 370 C (IF NDIST=1); NOT USED IF NDIST=2 MAN 380 C MAN 390 C ** NOTE ** A CONSTANT COLUMN OF 1@S CAN BE GENERATED MAN 400 C USING CLMEAN = CLMVAR = 1.0 AND NDIST = 1, MAN 410 C FOR EXAMPLE. MAN 420 C MAN 430 C NREPS -- NUMBER OF TIMES EACH UNIQUE ROW IS TO MAN 440 C APPEAR IN THE GENERATED DESIGN MATRIX MAN 450 C NDEG -- NUMBER OF ZERO RESIDUALS OUTSIDE OF THE MAN 460 C BASIS MAN 470 C RLOSS -- NUMBER OF DEPENDENT COLUMNS ADDED MAN 480 C NIND -- VALUE SUCH THAT THE DEPENDENT COLUMNS ARE MAN 490 C GENERATED BY SUMMING NIND SUCCESSIVE MAN 500 C COLUMNS MAN 510 C YDIST -- INTEGER WHICH SPECIFIES THE DISTRIBUTION MAN 520 C OF (NONZERO) RESIDUALS ASSOCIATED WITH MAN 530 C THE NONACTIVE CONSTRAINTS MAN 540 C = 0 IF NORMAL MAN 550 C = 1 IF UNIFORM MAN 560 C YMEAN -- MEAN OF RESIDUAL DISTRIBUTION (IF YDIST=0), MAN 570 C OR LOWER LIMIT (IF YDIST=1) MAN 580 C YVAR -- STANDARD DEVIATION OF RESIDUAL DISTRIBUTION MAN 590 C (IF YDIST=0), OR UPPER LIMIT (IF YDIST=1) MAN 600 C ISEED -- THE RANDOM NUMBER SEED USED TO MAN 610 C INITIATE THE RANDOM NUMBER GENERATOR MAN 620 C OUTPUT -- OUTPUT DEVISE UNIT TO WHICH THE GENERATED PROBLEM MAN 630 C WILL BE WRITTEN. MAN 640 C ** NOTE ** IF NXGEN,NDIST OR YDIST ARE OUTSIDE THE MAN 650 C THE INDICATED RANGES, THE REQUIRED DISTRIBUTION MAN 660 C IS DEFAULTED TO THAT OF A NORMAL. MAN 670 C IF ISEED IS SET TO ZERO BY THE USER MAN 680 C IT IS THEN DEFAULTED TO ONE ACCEPTABLE TO MAN 690 C THE RANDOM NUMBER GENERATOR. MAN 700 C MAN 710 C MAN 720 C OUTPUT OF THIS DRIVER' MAN 730 C ONE LINE WITH FORMAT (2I5) WHICH SPECIFIES NOBS AND NPAR, THEN MAN 740 C AS MANY LINES AS IS NECESSARY TO PRINT MAN 750 C THE AUGMENTED MATRIX CONTAINING IN ITS FIRST COLUMN THE RIGHT MAN 760 C HAND SIDE AND IN THE NEXT NPAR COLUMNS THE ELEMENTS OF THE MAN 770 C A-MATRIX FOR ALL NOBS ROWS. THIS MATRIX IS WRITTEN WITH MAN 780 C FORMAT (5E14.8). MAN 790 C MAN 800 C MAN 810 REAL A(25,400), RHS(400), CLMEAN(25), CLMVAR(25), XOPT(25) MAN 820 INTEGER NDIST(25), IACT(25), YDIST, RLOSS, OUTPUT MAN 830 C MAN 840 C READ (FROM UNIT 5) THE DATA NECESSARY TO GENERATE THE PROBLEM MAN 850 C MAN 860 READ (5,99999) NOBS, NPAR, NXGEN, NREPS, NDEG, RLOSS, NIND, OUTPUTMAN 870 NBAS = NPAR - RLOSS MAN 880 IF (NXGEN.NE.0) GO TO 10 MAN 890 READ (5,99998) (XOPT(I),I=1,NBAS) MAN 900 GO TO 20 MAN 910 10 READ (5,99998) XMEAN, XVAR MAN 920 20 READ (5,99997) (NDIST(I),I=1,NBAS) MAN 930 READ (5,99998) (CLMEAN(I),I=1,NBAS) MAN 940 READ (5,99998) (CLMVAR(I),I=1,NBAS) MAN 950 READ (5,99996) YDIST, YMEAN, YVAR MAN 960 READ (5,99995) ISEED MAN 970 IF (ISEED.NE.0) GO TO 30 MAN 980 ISEED = 2147483646 MAN 990 C MAN 1000 C GENERATE THE PROBLEM MAN 1010 C MAN 1020 30 CALL L1GNR(NOBS, NPAR, NXGEN, XOPT, XMEAN, XVAR, NDIST, CLMEAN, MAN 1030 * CLMVAR, NREPS, NDEG, RLOSS, NIND, YDIST, YMEAN, YVAR, A, RHS, MAN 1040 * IACT, SUMRES, ISEED, IERR) MAN 1050 C MAN 1060 C CHECK TO SEE THAT THE PROBLEM WAS GENERATED CORRECTLY MAN 1070 C MAN 1080 IF (IERR.NE.0) GO TO 50 MAN 1090 C MAN 1100 C WRITE OUT PROBLEM ON UNIT DEVICE = OUTPUT MAN 1110 C MAN 1120 WRITE (OUTPUT,99994) NOBS, NPAR MAN 1130 DO 40 J=1,NOBS MAN 1140 WRITE (OUTPUT,99993) RHS(J), (A(I,J),I=1,NPAR) MAN 1150 40 CONTINUE MAN 1160 GO TO 60 MAN 1170 C MAN 1180 C WRITE OUT ERROR WARNING MAN 1190 C MAN 1200 50 WRITE (6,99992) MAN 1210 60 STOP MAN 1220 99999 FORMAT (8I5) MAN 1230 99998 FORMAT (8F10.0) MAN 1240 99997 FORMAT (25I5) MAN 1250 99996 FORMAT (I5, 2F10.0) MAN 1260 99995 FORMAT (1I10) MAN 1270 99994 FORMAT (2I5) MAN 1280 99993 FORMAT (5E14.8) MAN 1290 99992 FORMAT (54H PROBLEM NOT PRINTED DUE TO ERROR IN DATA SPECIFICATIO,MAN 1300 * 1HN) MAN 1310 END MAN 1320 SUBROUTINE L1GNR(NOBS, NPAR, NXGEN, XOPT, XMEAN, XVAR, NDIST, L1G 10 * CLMEAN, CLMVAR, NREPS, NDEG, RLOSS, NIND, YDIST, YMEAN, YVAR, A, * RHS, IACT, SUMRES, ISEED, IERR) C C C ------------------------------------------------------------- C C DESCRIPTION' C C THIS SUBROUTINE GENERATES DATA SETS WHICH CAN BE C USED FOR TESTING L1-APPROXIMATION (LEAST ABSOLUTE C DEVIATION CURVE-FITTING) COMPUTER CODES. THE USER C CAN SPECIFY C C * PROBLEM SIZE C * SOLUTION VECTOR C * STATISTICAL DISTRIBUTION FOR COLUMN ELEMENTS C * ROW REPETITIONS C * DEGENERACY C * RANK LOSS C * STATISTICAL DISTRIBUTION FOR RESIDUALS C C ON OUTPUT, X = XOPT IS AN OPTIMAL SOLUTION TO THE C GENERATED L1-APPROXIMATION PROBLEM' C C C MINIMIZE // X*A - RHS // , C C C WHERE //...// IS THE L1-NORM AND THE TRANSPOSE OF C A REPRESENTS THE DESIGN MATRIX. FURTHERMORE, C X = XOPT IS GUARANTEED TO BE THE UNIQUE OPTIMAL C SOLUTION IF RLOSS=0. C C INPUT ARGUMENTS' C C NOBS -- NUMBER OF OBSERVATIONS C NPAR -- NUMBER OF PARAMETERS C NXGEN -- INDICATES WHETHER SOLUTION X IS SPECIFIED C AS INPUT OR GENERATED C = 0 IF SPECIFIED AS INPUT C = 1 IF RANDOMLY GENERATED FROM NORMAL C = 2 IF RANDOMLY GENERATED FROM UNIFORM C XOPT -- REAL ARRAY OF DIMENSION NPAR - RLOSS WHICH C SPECIFIES THE OPTIMAL SOLUTION IF NXGEN IS C SET EQUAL TO ZERO C XMEAN -- MEAN OF NORMAL DISTRIBUTION IF NXGEN=1, C LOWER LIMIT OF UNIFORM DISTRIBUTION IF C NXGEN=2 C XVAR -- STANDARD DEVIATION OF NORMAL DISTRIBUTION C IF NXGEN=1, UPPER LIMIT OF UNIFORM DISTRI- C BUTION IF NXGEN=2 C NDIST -- ARRAY OF DIMENSION NPAR - RLOSS WHICH C SPECIFIES THE STATISTICAL DISTRIBUTION C FOR EACH COLUMN OF THE DESIGN MATRIX C = 0 IF NORMAL C = 1 IF UNIFORM C = 2 IF ZERO-ONE C CLMEAN -- ARRAY OF DIMENSION NPAR - RLOSS WHICH C SPECIFIES FOR EACH COLUMN THE MEAN (IF C NDIST=0), LOWER LIMIT (IF NDIST=1), OR C PROPORTION P OF ZERO ENTRIES (IF NDIST=2) C CLMVAR -- ARRAY OF DIMENSION NPAR - RLOSS WHICH C SPECIFIES FOR EACH COLUMN THE STANDARD C DEVIATION (IF NDIST=0), THE UPPER LIMIT C (IF NDIST=1); NOT USED IF NDIST=2 C C ** NOTE ** A CONSTANT COLUMN OF 1@S CAN BE GENERATED C USING CLMEAN = CLMVAR = 1.0 AND NDIST = 1, C FOR EXAMPLE. C C NREPS -- NUMBER OF TIMES EACH UNIQUE ROW IS TO C APPEAR IN THE GENERATED DESIGN MATRIX C NDEG -- NUMBER OF ZERO RESIDUALS OUTSIDE OF THE C BASIS C RLOSS -- NUMBER OF DEPENDENT COLUMNS ADDED C NIND -- VALUE SUCH THAT THE DEPENDENT COLUMNS ARE C GENERATED BY SUMMING NIND SUCCESSIVE C COLUMNS C YDIST -- INTEGER WHICH SPECIFIES THE DISTRIBUTION C OF (NONZERO) RESIDUALS ASSOCIATED WITH C THE NONACTIVE CONSTRAINTS C = 0 IF NORMAL C = 1 IF UNIFORM C YMEAN -- MEAN OF RESIDUAL DISTRIBUTION (IF YDIST=0), C OR LOWER LIMIT (IF YDIST=1) C YVAR -- STANDARD DEVIATION OF RESIDUAL DISTRIBUTION C (IF YDIST=0), OR UPPER LIMIT (IF YDIST=1) C ISEED -- THE RANDOM NUMBER SEED USED TO C INITIATE THE RANDOM NUMBER GENERATOR C C ** NOTE ** IF NXGEN, NDIST OR YDIST ARE OUTSIDE THE C INDICATED RANGES, THE REQUIRED DISTRIBUTION C IS DEFAULTED TO THAT OF A NORMAL. C C OUTPUT ARGUMENTS' C C XOPT -- OPTIMAL SOLUTION VECTOR OF DIMENSION NPAR C A -- GENERATED TRANSPOSE OF THE DESIGN MATRIX, WITH C NPAR ROWS AND NOBS COLUMNS C RHS -- RIGHT HAND SIDE (OR Y) VECTOR OF DIMENSION C NOBS C IACT -- VECTOR OF DIMENSION NPAR - RLOSS WHICH C CONTAINS INDICES OF ACTIVE CONSTRAINT ROWS C SUMRES -- OBJECTIVE FUNCTION VALUE, OR SUM OF C ABSOLUTE VALUES OF RESIDUALS C ISEED -- THE RANDOM NUMBER SEED AVAILABLE C UPON TERMINATION OF GENERATION PROCESS C IERR -- ERROR FLAG, UPON RETURN C = 0 NORMAL EXECUTION C = 1 FATAL ERROR C C RESTRICTIONS' C C 1 .LE. NOBS .LE. 400 C 1 .LE. NPAR .LE. 25 C 1 .LE. NREPS .LE. NOBS C 0 .LE. NDEG .LE. NOBS-NPAR+RLOSS-2 C 0 .LE. RLOSS .LE. NPAR-1 C 0 .LE. NIND .LE. NPAR-RLOSS C NOBS+RLOSS-NPAR-NDEG MUST BE EVEN C NREPS MUST DIVIDE NOBS C C SUBROUTINES CALLED' C C MOD -- FORTRAN LIBRARY ROUTINE C NORRAN -- PSEUDO-RANDOM NUMBER GENERATOR FOR NORMAL C DISTRIBUTION WITH MEAN=0, STANDARD DEVIATION=1 C RAND -- PSEUDO-RANDOM NUMBER GENERATOR FOR UNIFORM C DISTRIBUTION OVER (0,1) C BINRAN -- PSEUDO-RANDOM NUMBER GENERATOR FOR ZERO-ONE C VARIABLES C RANPER -- PRODUCES A RANDOM PERMUTATION OF INTEGERS C 1 THROUGH N C C LANGUAGE' ANSI FORTRAN C C REFERENCE' C C K.L.HOFFMAN : D.R.SHIER, A TEST PROBLEM GENERATOR FOR C DISCRETE LINEAR L1 APPROXIMATION PROBLEMS@ C C WRITTEN BY C C KARLA L. HOFFMAN C OPERATIONS RESEARCH DIVISION C NATIONAL BUREAU OF STANDARDS C WASHINGTON, DC 20234 C C NOVEMBER, 1978 C C ------------------------------------------------------------- C C DIMENSION XOPT(25), NDIST(25), CLMEAN(25), CLMVAR(25), A(25,400), * RHS(400), IACT(25) DIMENSION B(25,400), Y(400), WORK(400) INTEGER YDIST, RLOSS, HALF DOUBLE PRECISION PROD, BB, SUM, SUM1 C C DEFINE OUTPUT UNIT FOR ERROR MESSAGES C IOUT = 6 C C CHECK CONSISTENCY OF INPUT C IERR = 0 IF (NOBS.LT.1 .OR. NOBS.GT.400) GO TO 400 IF (NPAR.LT.1 .OR. NPAR.GT.25) GO TO 410 IF (NREPS.LT.1 .OR. NREPS.GT.NOBS) GO TO 420 IF (RLOSS.LT.0 .OR. RLOSS.GE.NPAR) GO TO 430 NBAS = NPAR - RLOSS IF (NIND.LT.0 .OR. NIND.GT.NBAS) GO TO 440 NONBAS = NOBS - NBAS - NDEG IF (NDEG.LT.0 .OR. NDEG.GT.(NONBAS-2)) GO TO 450 IF (MOD(NONBAS,2).NE.0) GO TO 460 IF (MOD(NOBS,NREPS).NE.0) GO TO 470 C C GENERATE X-VECTOR IF NOT SPECIFIED C NBAS1 = NBAS + 1 IF (RLOSS.EQ.0) GO TO 20 DO 10 I=NBAS1,NPAR XOPT(I) = 0.0 10 CONTINUE 20 IF (NXGEN.EQ.0) GO TO 60 IF (NXGEN.EQ.2) GO TO 40 CALL NORRAN(NBAS, XOPT, ISEED) DO 30 I=1,NBAS XOPT(I) = XOPT(I)*XVAR + XMEAN 30 CONTINUE GO TO 60 40 CALL UNIRAN(NBAS, XOPT, ISEED) DO 50 I=1,NBAS XOPT(I) = XOPT(I)*(XVAR-XMEAN) + XMEAN 50 CONTINUE C C GENERATE THE DESIGN MATRIX OF DIFFERENT OBSERVATIONS (NO C REPETITIONS) C 60 NDIFOB = NOBS/NREPS DO 120 I=1,NBAS CLMV = CLMVAR(I) CLMM = CLMEAN(I) IF (NDIST(I).EQ.1) GO TO 80 IF (NDIST(I).EQ.2) GO TO 100 CALL NORRAN(NDIFOB, WORK, ISEED) DO 70 J=1,NDIFOB A(I,J) = WORK(J)*CLMV + CLMM 70 CONTINUE GO TO 120 80 CALL UNIRAN(NDIFOB, WORK, ISEED) DO 90 J=1,NDIFOB A(I,J) = WORK(J)*(CLMV-CLMM) + CLMM 90 CONTINUE GO TO 120 100 CALL BINRAN(NDIFOB, CLMM, 1, WORK, ISEED) DO 110 J=1,NDIFOB A(I,J) = WORK(J) 110 CONTINUE 120 CONTINUE C C CREATE THE REPETITIONS C IF (NREPS.EQ.1) GO TO 160 NREPS1 = NREPS - 1 DO 150 J=1,NDIFOB JJ = J DO 140 K=1,NREPS1 JJ = JJ + NDIFOB DO 130 I=1,NBAS A(I,JJ) = A(I,J) 130 CONTINUE 140 CONTINUE 150 CONTINUE C C PERMUTE ALL BUT THE FIRST NBAS ROWS OF THE DESIGN MATRIX C 160 K = NOBS - NBAS CALL RANPER(K, WORK, ISEED) DO 190 J=1,NOBS IF (J.LE.NBAS) K = J IF (J.LE.NBAS) GO TO 170 L = J - NBAS K = IFIX(WORK(L)) + NBAS 170 DO 180 I=1,NBAS B(I,K) = A(I,J) 180 CONTINUE 190 CONTINUE C C ADD THE DEPENDENT COLUMNS C IF (RLOSS.EQ.0) GO TO 230 L = 1 DO 220 I=NBAS1,NPAR L1 = L + NIND - 1 DO 210 J=1,NOBS BB = 0.0D0 DO 200 K=L,L1 BB = BB + DBLE(B(K,J)) 200 CONTINUE B(I,J) = BB 210 CONTINUE L = L + 1 220 CONTINUE C C CREATE FILL-IN ROW TO ENSURE OPTIMALITY OF BASIS C 230 HALF = NONBAS/2 L = NBAS1 L1 = NBAS + HALF L2 = L1 + 1 L4 = L1 + HALF L3 = L4 - 1 DO 260 I=1,NPAR SUM = 0.0D0 SUM1 = 0.0D0 DO 240 J=L,L1 SUM = SUM + DBLE(B(I,J)) 240 CONTINUE DO 250 J=L2,L3 SUM1 = SUM1 + DBLE(B(I,J)) 250 CONTINUE B(I,L4) = SUM - SUM1 260 CONTINUE C C CALCULATE THE RIGHT HAND SIDE BEFORE RESIDUALS ADDED C DO 280 J=1,NOBS PROD = 0.0D0 DO 270 I=1,NBAS PROD = PROD + (DBLE(B(I,J))*DBLE(XOPT(I))) 270 CONTINUE Y(J) = PROD 280 CONTINUE C C CALCULATE THE RESIDUALS C IF (YDIST.EQ.1) GO TO 300 CALL NORRAN(NONBAS, WORK, ISEED) DO 290 J=1,NONBAS WORK(J) = WORK(J)*YVAR + YMEAN 290 CONTINUE GO TO 320 300 CALL UNIRAN(NONBAS, WORK, ISEED) DO 310 J=1,NONBAS WORK(J) = WORK(J)*(YVAR-YMEAN) + YMEAN 310 CONTINUE C C FIND SUM OF RESIDUALS AND ADD RESIDUALS TO RIGHT HAND SIDE C 320 L = NBAS1 L1 = NOBS - NDEG HALF = L + HALF PROD = 0.0D0 DO 330 I=1,NONBAS PROD = PROD + DBLE(ABS(WORK(I))) 330 CONTINUE SUMRES = PROD K = 0 DO 350 J=L,L1 K = K + 1 IF (J.GE.HALF) GO TO 340 Y(J) = Y(J) + DBLE(ABS(WORK(K))) GO TO 350 340 Y(J) = Y(J) - DBLE(ABS(WORK(K))) 350 CONTINUE C C PERMUTE ALL ROWS OF THE DESIGN MATRIX C CALL RANPER(NOBS, WORK, ISEED) DO 370 J=1,NOBS K = WORK(J) RHS(K) = Y(J) IF (J.LE.NBAS) IACT(J) = K DO 360 I=1,NPAR A(I,K) = B(I,J) 360 CONTINUE 370 CONTINUE C C SORT ACTIVE CONSTRAINTS INTO ASCENDING ORDER C LBAS = NBAS - 1 IF (LBAS.EQ.0) GO TO 490 DO 390 J=1,LBAS JJ = J + 1 DO 380 I=JJ,NBAS IF (IACT(J).LE.IACT(I)) GO TO 380 ITEMP = IACT(J) IACT(J) = IACT(I) IACT(I) = ITEMP 380 CONTINUE 390 CONTINUE GO TO 490 400 WRITE (IOUT,99999) GO TO 480 410 WRITE (IOUT,99998) GO TO 480 420 WRITE (IOUT,99997) GO TO 480 430 WRITE (IOUT,99996) GO TO 480 440 WRITE (IOUT,99995) GO TO 480 450 WRITE (IOUT,99994) GO TO 480 460 WRITE (IOUT,99993) GO TO 480 470 WRITE (IOUT,99992) 480 IERR = 1 490 RETURN 99999 FORMAT (37H FATAL ERROR *** NOBS IS OUT OF RANGE) 99998 FORMAT (37H FATAL ERROR *** NPAR IS OUT OF RANGE) 99997 FORMAT (38H FATAL ERROR *** NREPS IS OUT OF RANGE) 99996 FORMAT (38H FATAL ERROR *** RLOSS IS OUT OF RANGE) 99995 FORMAT (37H FATAL ERROR *** NIND IS OUT OF RANGE) 99994 FORMAT (37H FATAL ERROR *** NDEG IS OUT OF RANGE) 99993 FORMAT (49H FATAL ERROR *** NOBS+RLOSS-NPAR-NDEG IS NOT EVEN) 99992 FORMAT (43H FATAL ERROR *** NREPS DOES NOT DIVIDE NOBS) END SUBROUTINE UNIRAN(N, RAND, IX) UNI 10 C PORTABLE RANDOM NUMBER GENERATOR C USING THE RECURSION C IX = IX*A MOD P C THIS FORTRAN IMPLEMENTATION OF A RANDOM GENERATOR C COMES FROM LINUS SCHRAGE AND IS PUBLISHED IN C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, JUNE, 1979. C C SOME COMPILERS, E.G., THE HP3000, REQUIRE THE FOLLOWING C DECLARATION TO BE INTEGER*4 INTEGER A, P, IX, B15, B16, XHI, XALO, LEFTLO, FHI, K DIMENSION RAND(1) C C C 7**5, 2**15, 2**16, 2**31-1 DATA A /16807/, B15 /32768/, B16 /65536/, P /2147483647/ C DO 10 I=1,N C GET 15 HI ORDER BITS OF IX XHI = IX/B16 C GET 16 LO BITS OF IX AND FORM LO PRODUCT XALO = (IX-XHI*B16)*A C GET 15 HI ORDER BITS OF LO PRODUCT LEFTLO = XALO/B16 C FORM THE 31 HIGHEST BITS OF FULL PRODUCT FHI = XHI*A + LEFTLO C GET OVERFLO PAST 31ST BIT OF FULL PRODUCT K = FHI/B15 C ASSEMBLE ALL THE PARTS AND PRESUBTRACT P C THE PARENTHESES ARE ESSENTIAL IX = (((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16) + K C ADD P BACK IN IF NECESSARY IF (IX.LT.0) IX = IX + P C MULTIPLY BY 1/(2**31-1) RAND(I) = FLOAT(IX)*4.656612875E-10 10 CONTINUE RETURN END SUBROUTINE NORRAN(N, X, ISEED) NOR 10 C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE THE NORMAL (GAUSSIAN) C DISTRIBUTION WITH MEAN = 0 AND STANDARD DEVIATION = 1. C THIS DISTRIBUTION IS DEFINED FOR ALL X AND HAS C THE PROBABILITY DENSITY FUNCTION C F(X) = (1/SQRT(2*PI))*EXP(-X*X/2). C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE NORMAL DISTRIBUTION C WITH MEAN = 0 AND STANDARD DEVIATION = 1. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--ALOG, SQRT, SIN, COS. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C METHOD--BOX-MULLER ALGORITHM. C REFERENCES--BOX AND MULLER, A NOTE ON THE GENERATION C OF RANDOM NORMAL DEVIATES@, JOURNAL OF THE C ASSOCIATION FOR COMPUTING MACHINERY, 1958, C PAGES 610-611. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 33-34. C --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, C 1964, PAGE 39. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 40-111. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE 301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --JULY 1976. C C --------------------------------------------------------------------- C DIMENSION X(1) DIMENSION Y(2) DATA PI /3.14159265358979/ C IPR = 6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF (N.LT.1) GO TO 10 GO TO 20 10 WRITE (IPR,99999) WRITE (IPR,99998) N RETURN 20 CONTINUE C C -----START POINT----------------------------------------------------- C C GENERATE N UNIFORM (0,1) RANDOM NUMBERS; C THEN GENERATE 2 ADDITIONAL UNIFORM (0,1) RANDOM NUMBERS C (TO BE USED BELOW IN FORMING THE N-TH NORMAL C RANDOM NUMBER WHEN THE DESIRED SAMPLE SIZE N C HAPPENS TO BE ODD). C CALL UNIRAN(N, X, ISEED) CALL UNIRAN(2, Y, ISEED) C C GENERATE N NORMAL RANDOM NUMBERS C USING THE BOX-MULLER METHOD. C DO 50 I=1,N,2 IP1 = I + 1 U1 = X(I) IF (I.EQ.N) GO TO 30 U2 = X(IP1) GO TO 40 30 U2 = Y(2) 40 ARG1 = -2.0*ALOG(U1) ARG2 = 2.0*PI*U2 SQRT1 = SQRT(ARG1) Z1 = SQRT1*COS(ARG2) Z2 = SQRT1*SIN(ARG2) X(I) = Z1 IF (I.EQ.N) GO TO 50 X(IP1) = Z2 50 CONTINUE C RETURN 99999 FORMAT (1H , 49H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO T, * 42HHE NORRAN SUBROUTINE IS NON-POSITIVE *****) 99998 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , I8, 6H *****) END SUBROUTINE BINRAN(N, P, NPAR, X, ISEED) BIN 10 C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION BERNOULLI PROBABILITY C PARAMETER = P, C AND INTEGER @NUMBER OF BERNOULLI TRIALS@ C PARAMETER = NPAR. C THE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = NPAR*P C AND STANDARD DEVIATION = SQRT(NPAR*P*(1-P)). C THIS DISTRIBUTION IS DEFINED FOR ALL C DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY) C AND NPAR (INCLUSIVELY). C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = C(NPAR,X) * P**X * (1-P)**(NPAR-X). C WHERE C(NPAR,X) IS THE COMBINATORIAL FUNCTION C EQUALING THE NUMBER OF COMBINATIONS OF NPAR ITEMS C TAKEN X AT A TIME. C THE BINOMIAL DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF C SUCCESSES IN NPAR BERNOULLI (0,1) C TRIALS WHERE THE PROBABILITY OF SUCCESS C IN A SINGLE TRIAL = P. C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C --P = THE SINGLE PRECISION VALUE C OF THE BERNOULLI PROBABILITY C PARAMETER FOR THE BINOMIAL C DISTRIBUTION. C P SHOULD BE BETWEEN C 0.0 (EXCLUSIVELY) AND C 1.0 (EXCLUSIVELY). C --NPAR = THE INTEGER VALUE C OF THE @NUMBER OF BERNOULLI TRIALS@ C PARAMETER. C NPAR SHOULD BE A POSITIVE INTEGER. C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE BINOMIAL DISTRIBUTION C WITH BERNOULLI PROBABILITY PARAMETER = P C AND NUMBER OF BERNOULLI TRIALS PARAMETER = NPAR. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C --NPAR SHOULD BE A POSITIVE INTEGER. C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN, GEORAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT C FROM THIS DISCRETE RANDOM NUMBER C GENERATOR MUST NECESSARILY BE A C SEQUENCE OF ***INTEGER*** VALUES, C THE OUTPUT VECTOR X IS SINGLE C PRECISION IN MODE. C X HAS BEEN SPECIFIED AS SINGLE C PRECISION SO AS TO CONFORM WITH THE DATAPAC C CONVENTION THAT ALL OUTPUT VECTORS FROM ALL C DATAPAC SUBROUTINES ARE SINGLE PRECISION. C THIS CONVENTION IS BASED ON THE BELIEF THAT C 1) A MIXTURE OF MODES (FLOATING POINT C VERSUS INTEGER) IS INCONSISTENT AND C AN UNNECESSARY COMPLICATION C IN A DATA ANALYSIS; AND C 2) FLOATING POINT MACHINE ARITHMETIC C (AS OPPOSED TO INTEGER ARITHMETIC) C IS THE MORE NATURAL MODE FOR DOING C DATA ANALYSIS. C REFERENCES--JOHNSON AND KOTZ, DISCRETE C DISTRIBUTIONS, 1969, PAGES 50-86. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGE 41. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 135-142. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 929. C --KENDALL AND STUART, THE ADVANCED THEORY OF C STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 120-125. C --MOOD AND GRABLE, INTRODUCTION TO THE THEORY C OF STATISTICS, EDITION 2, 1963, PAGES 64-69. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 39-40. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE' 301-921-2315 C ORIGINAL VERSION--NOVEMBER 1975. C C --------------------------------------------------------------------- C DIMENSION X(1) C IPR = 6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF (N.LT.1) GO TO 10 IF (P.LE.0.0 .OR. P.GE.1.0) GO TO 20 IF (NPAR.LT.1) GO TO 30 GO TO 40 10 WRITE (IPR,99999) WRITE (IPR,99995) N RETURN 20 WRITE (IPR,99998) WRITE (IPR,99996) P RETURN 30 WRITE (IPR,99997) WRITE (IPR,99995) NPAR RETURN 40 CONTINUE C C -----START POINT----------------------------------------------------- C C CHECK ON THE MAGNITUDE OF P, C AND BRANCH TO THE FASTER C GENERATION METHOD ACCORDINGLY. C IF (P.LT.0.1) GO TO 70 C C IF P IS MODERATE OR LARGE, C GENERATE N BINOMIAL RANDOM NUMBERS C USING THE REJECTION METHOD. C DO 60 I=1,N ISUM = 0 DO 50 J=1,NPAR CALL UNIRAN(1, U, ISEED) IF (U.LE.P) ISUM = ISUM + 1 50 CONTINUE X(I) = ISUM 60 CONTINUE RETURN C C IF P IS SMALL, C GENERATE N BINOMIAL NUMBERS C USING THE FACT THAT THE C WAITING TIME FOR 1 SUCCESS IN C BERNOULLI TRIALS HAS A C GEOMETRIC DISTRIBUTION. C 70 DO 100 I=1,N ISUM = 0 J = 1 80 CALL GEORAN(1, P, G, ISEED) IG = G + 0.5 ISUM = ISUM + IG + 1 IF (ISUM.GT.NPAR) GO TO 90 J = J + 1 GO TO 80 90 X(I) = J - 1 100 CONTINUE RETURN C 99999 FORMAT (1H , 49H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO T, * 42HHE BINRAN SUBROUTINE IS NON-POSITIVE *****) 99998 FORMAT (1H , 49H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO T, * 61HHE BINRAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL , * 5H*****) 99997 FORMAT (1H , 49H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO T, * 42HHE BINRAN SUBROUTINE IS NON-POSITIVE *****) 99996 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , E15.8, * 6H *****) 99995 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , I8, 6H *****) END SUBROUTINE GEORAN(N, P, X, ISEED) GEO 10 C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE GEOMETRIC DISTRIBUTION C WITH SINGLE PRECISION ^ERNOULLI PROBABILITY@ G C PARAMETER = P. C THE GEOMETRIC DISTRIBUTION USED C HEREIN HAS MEAN = (1-P)/P C AND STANDARD DEVIATION = SQRT((1-P)/(P*P))). C THIS DISTRIBUTION IS DEFINED FOR C ALL NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... . C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = P * (1-P)**X. C THE GEOMETRIC DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF FAILURES C BEFORE OBTAINING 1 SUCCESS IN AN C INDEFINITE SEQUENCE OF BERNOULLI (0,1) C TRIALS WHERE THE PROBABILITY OF SUCCESS C IN A SINGLE TRIAL = P. C INPUT ARGUMENTS--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. C --P = THE SINGLE PRECISION VALUE C OF THE ^ERNOULLI PROBABILITY@ G C PARAMETER FOR THE GEOMETRIC C DISTRIBUTION. C P SHOULD BE BETWEEN C 0.0 (EXCLUSIVELY) AND C 1.0 (EXCLUSIVELY). C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE GEOMETRIC DISTRIBUTION C WITH ^ERNOULLI PROBABILITY@ PARAMETER = P. G C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT C FROM THIS DISCRETE RANDOM NUMBER C GENERATOR MUST NECESSARILY BE A C SEQUENCE OF ***INTEGER*** VALUES, C THE OUTPUT VECTOR X IS SINGLE C PRECISION IN MODE. C X HAS BEEN SPECIFIED AS SINGLE C PRECISION SO AS TO CONFORM WITH THE DATAPAC C CONVENTION THAT ALL OUTPUT VECTORS FROM ALL C DATAPAC SUBROUTINES ARE SINGLE PRECISION. C THIS CONVENTION IS BASED ON THE BELIEF THAT C 1) A MIXTURE OF MODES (FLOATING POINT C VERSUS INTEGER) IS INCONSISTENT AND C AN UNNECESSARY COMPLICATION C IN A DATA ANALYSIS; AND C 2) FLOATING POINT MACHINE ARITHMETIC C (AS OPPOSED TO INTEGER ARITHMETIC) C IS THE MORE NATURAL MODE FOR DOING C DATA ANALYSIS. C REFERENCES--TOCHER, THE ART OF SIMULATION, C 1963, PAGES 14-15. C --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, C 1964, PAGE 36. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 155-157, 210. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 929. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE' 301-921-2315 C ORIGINAL VERSION--NOVEMBER 1975. C C --------------------------------------------------------------------- C DIMENSION X(1) C IPR = 6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF (N.LT.1) GO TO 10 IF (P.LE.0.0 .OR. P.GE.1.0) GO TO 20 GO TO 30 10 WRITE (IPR,99999) WRITE (IPR,99996) N RETURN 20 WRITE (IPR,99998) WRITE (IPR,99997) P RETURN 30 CONTINUE C C -----START POINT----------------------------------------------------- C C GENERATE N UNIFORM (0,1) RANDOM NUMBERS; C CALL UNIRAN(N, X, ISEED) C C GENERATE N GEOMETRIC RANDOM NUMBERS C USING THE PERCENT POINT FUNCTION TRANSFORMATION METHOD. C DO 40 I=1,N IF (X(I).LE.0.0 .OR. X(I).GE.1.0) GO TO 40 ARG1 = 1.0 - X(I) ARG2 = 1.0 - P ANUM = ALOG(ARG1) ADEN = ALOG(ARG2) RATIO = ANUM/ADEN IRATIO = RATIO X(I) = IRATIO ARATIO = IRATIO IF (ARATIO.EQ.RATIO) X(I) = IRATIO - 1 40 CONTINUE C RETURN 99999 FORMAT (1H , 49H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO T, * 42HHE GEORAN SUBROUTINE IS NON-POSITIVE *****) 99998 FORMAT (1H , 49H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO T, * 61HHE GEORAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL , * 5H*****) 99997 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , E15.8, * 6H *****) 99996 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , I8, 6H *****) END SUBROUTINE RANPER(N, X, ISEED) RAN 10 C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM PERMUTATION OF SIZE N C OF THE VALUES 1.0, 2.0, 3.0, ..., N-1, N. C INPUT ARGUMENTS--N = THE DESIRED INTEGER SIZE C OF THE RANDOM 1 TO N PERMUTATION. C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM PERMUTATION WILL BE PLACED. C OUTPUT--A RANDOM PERMUTATION OF SIZE N C OF THE VALUES 1.0, 2.0, 3.0, ..., N-1, N. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--ALGORITHM SUGGESTED BY DAN LOZIER, C NATIONAL BUREAU OF STANDARDS (205.01). C REFERENCES--NONE. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE' 301-921-2315 C ORIGINAL VERSION--JUNE 1972. C UPDATED --MAY 1974. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C C --------------------------------------------------------------------- C DIMENSION X(1) DIMENSION U(1) C AN = N IPR = 6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF (N.LT.1) GO TO 10 IF (N.EQ.1) GO TO 20 GO TO 30 10 WRITE (IPR,99999) WRITE (IPR,99997) N RETURN 20 WRITE (IPR,99998) X(1) = 1 RETURN 30 CONTINUE C C -----START POINT----------------------------------------------------- C DO 40 I=1,N X(I) = I 40 CONTINUE C DO 50 I=1,N CALL UNIRAN(1, U(1), ISEED) ADD = AN*U(1) + 1.0 IADD = ADD IF (IADD.LT.1) IADD = 1 IF (IADD.GT.N) IADD = N J = I + IADD IF (J.GT.N) J = J - N HOLD = X(J) X(J) = X(I) X(I) = HOLD 50 CONTINUE RETURN 99999 FORMAT (1H , 49H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO T, * 42HHE RANPER SUBROUTINE IS NON-POSITIVE *****) 99998 FORMAT (1H , 49H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGU, * 51HMENT TO THE RANPER SUBROUTINE HAS THE VALUE 1 *****) 99997 FORMAT (1H , 35H***** THE VALUE OF THE ARGUMENT IS , I8, 6H *****) END