*DECK UERTST C IMSL ROUTINE NAME - UERTST C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION C C USAGE - CALL UERTST (IER,NAME) C C ARGUMENTS IER - ERROR PARAMETER. (INPUT) C IER = I+J WHERE C I = 128 IMPLIES TERMINAL ERROR, C I = 64 IMPLIES WARNING WITH FIX, AND C I = 32 IMPLIES WARNING. C J = ERROR CODE RELEVANT TO CALLING C ROUTINE. C NAME - A SIX CHARACTER LITERAL STRING GIVING THE C NAME OF THE CALLING ROUTINE. (INPUT) C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN C ONTO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT C NUMBER CAN BE DETERMINED BY CALLING UGETIO AS C FOLLOWS.. CALL UGETIO(1,NIN,NOUT). C THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING C UGETIO AS FOLLOWS.. C NIN = 0 C NOUT = NEW OUTPUT UNIT NUMBER C CALL UGETIO(3,NIN,NOUT) C SEE THE UGETIO DOCUMENT FOR MORE DETAILS. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE UERTST (IER,NAME) C SPECIFICATIONS FOR ARGUMENTS INTEGER IER CHARACTER NAME*6 C SPECIFICATIONS FOR LOCAL VARIABLES CHARACTER*6 NAMSET,NAMEQ,IEQ*1 DATA NAMSET/'UERSET'/ DATA NAMEQ/' '/ C FIRST EXECUTABLE STATEMENT DATA LEVEL/4/,IEQDF/0/,IEQ/'='/ IF (IER.GT.999) GO TO 25 IF (IER.LT.-32) GO TO 55 IF (IER.LE.128) GO TO 5 IF (LEVEL.LT.1) GO TO 30 C PRINT TERMINAL MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,35) IER,NAME GO TO 30 5 IF (IER.LE.64) GO TO 10 IF (LEVEL.LT.2) GO TO 30 C PRINT WARNING WITH FIX MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAME GO TO 30 10 IF (IER.LE.32) GO TO 15 C PRINT WARNING MESSAGE IF (LEVEL.LT.3) GO TO 30 CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,45) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAME GO TO 30 15 CONTINUE C CHECK FOR UERSET CALL IF (NAME.NE.NAMSET) GO TO 25 LEVOLD = LEVEL LEVEL = IER IER = LEVOLD IF (LEVEL.LT.0) LEVEL = 4 IF (LEVEL.GT.4) LEVEL = 4 GO TO 30 25 CONTINUE IF (LEVEL.LT.4) GO TO 30 C PRINT NON-DEFINED MESSAGE CALL UGETIO(1,NIN,IOUNIT) IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER,NAMEQ,IEQ,NAME IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER,NAME 30 IEQDF = 0 RETURN 35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 40 FORMAT(36H *** WARNING WITH FIX ERROR (IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) 50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5, 1 20H) FROM IMSL ROUTINE ,A6,A1,A6) C SAVE P FOR P = R CASE C P IS THE PAGE NAME C R IS THE ROUTINE NAME 55 IEQDF = 1 NAMEQ = NAME 65 RETURN END *DECK UGETIO C IMSL ROUTINE NAME - UGETIO C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - TO RETRIEVE CURRENT VALUES AND TO SET NEW C VALUES FOR INPUT AND OUTPUT UNIT C IDENTIFIERS. C C USAGE - CALL UGETIO(IOPT,NIN,NOUT) C C ARGUMENTS IOPT - OPTION PARAMETER. (INPUT) C IF IOPT=1, THE CURRENT INPUT AND OUTPUT C UNIT IDENTIFIER VALUES ARE RETURNED IN NIN C AND NOUT, RESPECTIVELY. C IF IOPT=2, THE INTERNAL VALUE OF NIN IS C RESET FOR SUBSEQUENT USE. C IF IOPT=3, THE INTERNAL VALUE OF NOUT IS C RESET FOR SUBSEQUENT USE. C NIN - INPUT UNIT IDENTIFIER. C OUTPUT IF IOPT=1, INPUT IF IOPT=2. C NOUT - OUTPUT UNIT IDENTIFIER. C OUTPUT IF IOPT=1, INPUT IF IOPT=3. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - NONE REQUIRED C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C REMARKS EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT C OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT C IDENTIFIER VALUES. IF UGETIO IS CALLED WITH IOPT=2 OR C IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED. C SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE UGETIO(IOPT,NIN,NOUT) C SPECIFICATIONS FOR ARGUMENTS INTEGER IOPT,NIN,NOUT C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER NIND,NOUTD C FIRST EXECUTABLE STATEMENT IF (IOPT.EQ.3) GO TO 10 IF (IOPT.EQ.2) GO TO 5 IF (IOPT.NE.1) GO TO 9005 NIN = I1MACH(1) NOUT = I1MACH(2) GO TO 9005 5 NIND = NIN GO TO 9005 10 NOUTD = NOUT 9005 RETURN END *DECK MERFI C IMSL ROUTINE NAME - MERFI C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - INVERSE ERROR FUNCTION C C USAGE - CALL MERFI (P,Y,IER) C C ARGUMENTS P - INPUT VALUE IN THE EXCLUSIVE RANGE (-1.0,1.0) C Y - OUTPUT VALUE OF THE INVERSE ERROR FUNCTION C IER - ERROR PARAMETER (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES P LIES OUTSIDE THE LEGAL C RANGE. PLUS OR MINUS MACHINE INFINITY IS C GIVEN AS THE RESULT (SIGN IS THE SIGN OF C THE FUNCTION VALUE OF THE NEAREST LEGAL C ARGUMENT). C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MERFI (P,Y,IER) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y INTEGER IER C SPECIFICATIONS FOR LOCAL VARIABLES REAL A(65),H1,H2,H3,H4,X,SIGMA,Z,W,X3,X4,X5, * X6,B INTEGER N,IPP,L,LB2 DATA A(1),A(2),A(3),A(4),A(5),A(6),A(7),A(8),A(9), * A(10),A(11),A(12),A(13),A(14),A(15),A(16), * A(17),A(18),A(19),A(20),A(21),A(22),A(23) * /.99288537662,.12046751614, * .16078199342E-01,.26867044372E-02, * .49963473024E-03,.98898218599E-04, * .20391812764E-04,.4327271618E-05, * .938081413E-06,.206734721E-06, * .46159699E-07,.10416680E-07, * .2371501E-08,.543928E-09, * .125549E-09,.29138E-10, * .6795E-11,.1591E-11, * .374E-12,.88E-13, * .21E-13,.5E-14, * .1E-14/ DATA A(24),A(25),A(26),A(27),A(28),A(29),A(30), * A(31),A(32),A(33),A(34),A(35),A(36),A(37), * A(38),A(39),A(40) * /.91215880342E00,-.16266281868E-01, * .43355647295E-03,.21443857007E-03, * .2625751076E-05,-.302109105E-05, * -.12406062E-07,.62406609E-07, * -.540125E-09,-.142328E-08, * .34384E-10,.33585E-10, * -.1458E-11,-.81E-12, * .53E-13,.2E-13, * -.2E-14/ DATA A(41),A(42),A(43),A(44),A(45),A(46),A(47), * A(48),A(49),A(50),A(51),A(52),A(53),A(54), * A(55),A(56),A(57),A(58),A(59),A(60),A(61), * A(62),A(63),A(64),A(65) * /.95667970902,-.023107004309, * -.43742360975E-02,-.57650342265E-03, * -.10961022307E-04,.25108547025E-04, * .10562336068E-04,.275441233E-05, * .432484498E-06,-.20530337E-07, * -.43891537E-07,-.1768401E-07, * -.3991289E-08,-.186932E-09, * .272923E-09,.132817E-09, * .31834E-10,.1670E-11, * -.2036E-11,-.965E-12, * -.22E-12,-.1E-13, * .13E-13,.6E-14, * .1E-14/ DATA H1,H2,H3,H4/-1.5488130424, * 2.5654901231,-.55945763133, * 2.2879157163/ C FIRST EXECUTABLE STATEMENT X = P IER = 0 SIGMA = SIGN(1.,X) IF (.NOT.(X.GT.-1..AND.X.LT.1.)) GO TO 35 Z = ABS(X) IF(Z.GT..8) GO TO 20 W = Z*Z/.32-1. N = 22 IPP = 1 L = 1 5 LB2 = 1 X3 = 1. X4 = W X6 = A(IPP) 10 X6 = X6 + A(IPP+LB2) * X4 X5 = X4 * W * 2.-X3 X3 = X4 X4 = X5 LB2 = LB2 + 1 IF (LB2 .LE. N) GO TO 10 GO TO (15,30),L 15 Y = Z * X6 * SIGMA GO TO 9005 20 B = SQRT(-ALOG(1.-Z*Z)) IF (Z .GT. .9975) GO TO 25 W = H1*B+H2 IPP = 24 L = 2 N = 16 GO TO 5 25 W = H3 * B + H4 IPP = 41 N = 24 L = 2 GO TO 5 30 Y = B * X6 * SIGMA GO TO 9005 35 Y = SIGMA*R1MACH(2) IER = 129 9000 CONTINUE CALL UERTST(IER,'MERFI ') 9005 RETURN END *DECK MDNRIS C IMSL ROUTINE NAME - MDNRIS C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1981 C C PURPOSE - INVERSE STANDARD NORMAL (GAUSSIAN) C PROBABILITY DISTRIBUTION FUNCTION C C USAGE - CALL MDNRIS (P,Y,IER) C C ARGUMENTS P - INPUT VALUE IN THE EXCLUSIVE RANGE (0.0,1.0) C Y - OUTPUT VALUE OF THE INVERSE NORMAL (0,1) C PROBABILITY DISTRIBUTION FUNCTION C IER - ERROR PARAMETER (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES P LIES OUTSIDE THE LEGAL C RANGE. PLUS OR MINUS MACHINE INFINITY IS C GIVEN AS THE RESULT (SIGN IS THE SIGN OF C THE FUNCTION VALUE OF THE NEAREST LEGAL C ARGUMENT). C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - MERFI,UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDNRIS (P,Y,IER) C SPECIFICATIONS FOR ARGUMENTS REAL P,Y INTEGER IER C SPECIFICATIONS FOR LOCAL VARIABLES REAL D(25),A,B,W,H3,H4,X3,X4,X5,X6 REAL SIGMA,SQRT2,X DATA SQRT2/1.4142135623731/ DATA D(1)/.95667970902049/ DATA D(2)/-.02310700430907/ DATA D(3)/-.00437423609751/ DATA D(4)/-.00057650342265/ DATA D(5)/-.00001096102231/ DATA D(6)/.00002510854703/ DATA D(7)/.00001056233607/ DATA D(8)/.00000275441233/ DATA D(9)/.00000043248450/ DATA D(10)/-.00000002053034/ DATA D(11)/-.00000004389154/ DATA D(12)/-.00000001768401/ DATA D(13)/-.00000000399129/ DATA D(14)/-.00000000018693/ DATA D(15)/.00000000027292/ DATA D(16)/.00000000013282/ DATA D(17)/.00000000003183/ DATA D(18)/.00000000000167/ DATA D(19)/-.00000000000204/ DATA D(20)/-.00000000000097/ DATA D(21)/-.00000000000022/ DATA D(22)/-.00000000000001/ DATA D(23)/.00000000000001/ DATA D(24)/.00000000000001/ DATA D(25)/.00000000000000/ DATA H3/-.55945763132983/ DATA H4/2.2879157162634/ C FIRST EXECUTABLE STATEMENT IER = 0 IF (P .GT. 0.0 .AND. P .LT. 1.0) GO TO 5 IER = 129 SIGMA = SIGN(1.0,P) Y = SIGMA * R1MACH(2) GO TO 9000 5 IF(P.LE.R1MACH(4)) GO TO 10 X = 1.0 -(P + P) CALL MERFI (X,Y,IER) Y = -SQRT2 * Y GO TO 9005 C P TOO SMALL, COMPUTE Y DIRECTLY 10 A = P+P B = SQRT(-ALOG(A+(A-A*A))) W = H3 * B + H4 X3 = 1.0 X4 = W X6 = D(1) DO 15 I=2,25 X6 = X6 + D(I) * X4 X5 = X4 * W * 2. - X3 X3 = X4 X4 = X5 15 CONTINUE Y = B*X6 Y = -SQRT2 * Y GO TO 9005 9000 CONTINUE CALL UERTST(IER,'MDNRIS') 9005 RETURN END *DECK MDCH C IMSL ROUTINE NAME - MDCH C C----------------------------------------------------------------------- C C COMPUTER - CDC/SINGLE C C LATEST REVISION - JUNE 1, 1980 C C PURPOSE - CHI-SQUARED PROBABILITY DISTRIBUTION FUNCTION C C USAGE - CALL MDCH (CS,DF,P,IER) C C ARGUMENTS CS - INPUT VALUE FOR WHICH THE PROBABILITY IS C COMPUTED. CS MUST BE GREATER THAN OR EQUAL C TO ZERO. C DF - INPUT NUMBER OF DEGREES OF FREEDOM OF THE C CHI-SQUARED DISTRIBUTION. DF MUST BE GREATER C THAN OR EQUAL TO .5 AND LESS THAN OR EQUAL C TO 200,000. C P - OUTPUT PROBABILITY THAT A RANDOM VARIABLE C WHICH FOLLOWS THE CHI-SQUARED DISTRIBUTION C WITH DF DEGREES OF FREEDOM IS LESS THAN OR C EQUAL TO CS. C IER - ERROR PARAMETER. (OUTPUT) C TERMINAL ERROR C IER = 129 INDICATES THAT CS OR DF WAS C SPECIFIED INCORRECTLY. C WARNING ERROR C IER = 34 INDICATES THAT THE NORMAL PDF C WOULD HAVE PRODUCED AN UNDERFLOW. C C PRECISION/HARDWARE - SINGLE/ALL C C REQD. IMSL ROUTINES - H32/MDNOR,MERRC=ERFC,MGAMAD=DGAMMA,UERTST, C UGETIO C - H36,H48,H60/MDNOR,MERRC=ERFC,MGAMA=GAMMA, C UERTST,UGETIO C C NOTATION - INFORMATION ON SPECIAL NOTATION AND C CONVENTIONS IS AVAILABLE IN THE MANUAL C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP C C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. C C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN C APPLIED TO THIS CODE. NO OTHER WARRANTY, C EXPRESSED OR IMPLIED, IS APPLICABLE. C C----------------------------------------------------------------------- C SUBROUTINE MDCH (CS,DF,P,IER) C SPECIFICATIONS FOR ARGUMENTS REAL CS,DF,P C SPECIFICATIONS FOR LOCAL VARIABLES REAL A,Z,DGAM,EPS,W,W1,B,Z1,ONE,HALF,X,C REAL GAMMA,THRD,PT2 DATA EPS/1.E-14/,HALF/5.E-1/ DATA THRTEN/13.E0/,ONE/1.E0/ DATA THRD/.33333333333333E0/ DATA PT2/.22222222222222E0/ FUNC(W,A,Z)=W* EXP(A*ALOG(Z)-Z) C FIRST EXECUTABLE STATEMENT C TEST FOR INVALID INPUT VALUES IF (DF .GE. .5 .AND. DF .LE. 2.E5 .AND. CS .GE. 0.0) GO TO 5 IER=129 P=-R1MACH(1) GO TO 9000 5 IER=0 C SET P=0. IF CS IS LESS THAN OR C EQUAL TO 10.**(-12) IF (CS .GT. 1.E-12) GO TO 15 10 P=0.0 GO TO 9005 15 IF(DF.LE.100.) GO TO 20 C USE NORMAL DISTRIBUTION APPROXIMATION C FOR LARGE DEGREES OF FREEDOM IF(CS.LT.2.0) GO TO 10 X=((CS/DF)**THRD-(ONE-PT2/DF))/SQRT(PT2/DF) IF (X .GT. 5.0) GO TO 50 IF (X .LT. -11.313708) GO TO 55 CALL MDNOR (X,P) GO TO 9005 C INITIALIZATION FOR CALCULATION USING C INCOMPLETE GAMMA FUNCTION 20 IF (CS .GT. 200.) GO TO 50 A=HALF*DF Z=HALF*CS C = A DGAM = GAMMA(C) W=AMAX1(HALF*A,THRTEN) IF (Z .GE. W) GO TO 35 IF (DF .GT. 25. .AND. CS .LT. 2.) GO TO 10 C CALCULATE USING EQUATION NO. 6.5.29 W=ONE/(DGAM*A) W1=W DO 25 I=1,50 B=I W1=W1*Z/(A+B) IF (W1 .LE. EPS*W) GO TO 30 W=W+W1 25 CONTINUE 30 P=FUNC(W,A,Z) GO TO 9005 C CALCULATE USING EQUATION NO. 6.5.32 35 Z1=ONE/Z B=A-ONE W1=B*Z1 W=ONE+W1 DO 40 I=2,50 B=B-ONE W1=W1*B*Z1 IF (W1 .LE. EPS*W) GO TO 45 W=W+W1 40 CONTINUE 45 W=Z1*FUNC(W,A,Z) P=ONE-W/DGAM GO TO 9005 50 P=1.0 GO TO 9005 C WARNING ERROR - UNDERFLOW WOULD HAVE C OCCURRED 55 P=0.0 IER=34 9000 CONTINUE CALL UERTST (IER,'MDCH ') 9005 RETURN END