annotate libcruft/ranlib/genf.f @ 5915:b2e1be30c8e9 ss-2-9-7

[project @ 2006-07-28 18:08:56 by jwe]
author jwe
date Fri, 28 Jul 2006 18:08:56 +0000
parents df7c57a6639d
children 0e71ead7359d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 REAL FUNCTION genf(dfn,dfd)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2 C**********************************************************************
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 C REAL FUNCTION GENF( DFN, DFD )
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 C GENerate random deviate from the F distribution
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 C Function
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 C Generates a random deviate from the F (variance ratio)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 C distribution with DFN degrees of freedom in the numerator
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 C and DFD degrees of freedom in the denominator.
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 C Arguments
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 C DFN --> Numerator degrees of freedom
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 C (Must be positive)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 C REAL DFN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 C DFD --> Denominator degrees of freedom
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 C (Must be positive)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 C REAL DFD
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 C Method
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
30 C Directly generates ratio of chisquare variates
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
31 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
32 C**********************************************************************
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
33 C .. Scalar Arguments ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
34 REAL dfd,dfn
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
35 C ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
36 C .. Local Scalars ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
37 REAL xden,xnum
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
38 C ..
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
39 C JJV changed this code to call sgamma directly
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
40 C .. External Functions ..
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
41 C REAL genchi
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
42 C EXTERNAL genchi
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
43 REAL sgamma
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
44 EXTERNAL sgamma
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
45 C ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
46 C .. Executable Statements ..
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
47 IF (.NOT. (dfn.LE.0.0.OR.dfd.LE.0.0)) GO TO 10
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
48 WRITE (*,*) 'Degrees of freedom nonpositive in GENF - abort!'
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
49 WRITE (*,*) 'DFN value: ',dfn,'DFD value: ',dfd
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
50 STOP 'Degrees of freedom nonpositive in GENF - abort!'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
51
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
52 10 xnum = 2.0*sgamma(dfn/2.0)/dfn
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 C GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
55 xden = 2.0*sgamma(dfd/2.0)/dfd
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
56 C JJV changed constant so that it will not underflow at compile time
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
57 C JJV while not slowing generator by using double precision or logs.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
58 C IF (.NOT. (xden.LE. (1.0E-38*xnum))) GO TO 20
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
59 IF (.NOT. (xden.LE. (1.0E-37*xnum))) GO TO 20
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
60 WRITE (*,*) ' GENF - generated numbers would cause overflow'
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 WRITE (*,*) ' Numerator ',xnum,' Denominator ',xden
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
62 C JJV next 2 lines changed to maintain truncation of large deviates.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
63 C WRITE (*,*) ' GENF returning 1.0E38'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
64 C genf = 1.0E38
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
65 WRITE (*,*) ' GENF returning 1.0E37'
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
66 genf = 1.0E37
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
67 GO TO 30
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
68
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
69 20 genf = xnum/xden
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
70 30 RETURN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
71
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
72 END