C ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488. c Driving program to read and c process the data file integer i, j, ncol, nrow double precision emin, expect, percnt, pre, prt, table(200) c external fexact c Read table dimensions 10 read (*,*) nrow, ncol, expect, percnt, emin c Terminate on 0, 0 if (nrow.eq.0 .and. ncol.eq.0) go to 9000 c Read and output TABLE write (*,99998) do 20 i=1, nrow read (*,*) (table(i+(j-1)*nrow),j=1,ncol) write (*,99997) (table(i+(j-1)*nrow),j=1,ncol) 99997 format (1x, 10f7.0) 20 continue c call fexact (nrow, ncol, table, nrow, expect, percnt, emin, prt, & pre) c write (*,99999) prt, pre go to 10 c 99998 format (/, 2x, 'The contingency table for this problem is:') 99999 format (2x, 'PRT =', f8.6, ' PRE = ', f8.6) 9000 stop end c----------------------------------------------------------------------- c Name: FEXACT c c Purpose: Computes Fisher's exact test probabilities and a hybrid c approximation to Fisher exact test probabilities for a c contingency table using the network algorithm. c c Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT, c EMIN, PRT, PRE) c c Arguments: c NROW - The number of rows in the table. (Input) c NCOL - The number of columns in the table. (Input) c TABLE - NROW by NCOL matrix containing the contingency table. c (Input) c LDTABL - Leading dimension of TABLE exactly as specified in the c dimension statement in the calling program. (Input) c EXPECT - Expected value used in the hybrid algorithm for c deciding when to use asymptotic theory probabilities. c (Input) c If EXPECT .LE. 0.0 then asymptotic theory probabilities c are not used and Fisher exact test probabilities are c computed. Otherwise, if PERCNT or more of the cells in c the remaining table have estimated expected values of c EXPECT or more, with no remaining cell having expected c value less than EMIN, then asymptotic chi-squared c probabilities are used. See the algorithm section of the c manual document for details. Use EXPECT = 5.0 to obtain c the 'Cochran' condition. c PERCNT - Percentage of remaining cells that must have estimated c expected values greater than EXPECT before asymptotic c probabilities can be used. (Input) c See argument EXPECT for details. Use PERCNT = 80.0 to c obtain the 'Cochran' condition. c EMIN - Minimum cell estimated expected value allowed for c asymptotic chi-squared probabilities to be used. (Input) c See argument EXPECT for details. Use EMIN = 1.0 to c obtain the 'Cochran' condition. c PRT - Probability of the observed table for fixed marginal c totals. (Output) c PRE - Table p-value. (Output) c PRE is the probability of a more extreme table, where c 'extreme' is in a probabilistic sense. c If EXPECT .LT. 0 then the Fisher exact probability c is returned. Otherwise, an approximation to the c Fisher exact probability is computed based upon c asymptotic chi-squared probabilities for ``large'' c table expected values. The user defines ``large'' c through the arguments EXPECT, PERCNT, and EMIN. c c Remarks: c 1. For many problems one megabyte or more of workspace can be c required. If the environment supports it, the user should begin c by increasing the workspace used to 200,000 units. c c 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used c by STP may be changed by changing the line MULT = 30 below to c another value. c c 3. FEXACT may be converted to single precision by setting IREAL = 3, c and converting all DOUBLE PRECISION specifications (except the c specifications for RWRK, IWRK, and DWRK) to REAL. This will c require changing the names and specifications of the intrinsic c functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the c machine specific constants will need to be changed, and the name c DWRK will need to be changed to RWRK in the call to F2XACT. c c 4. Machine specific constants are specified and documented in F2XACT. c A missing value code is specified in both FEXACT and F2XACT. c c 5. Although not a restriction, is is not generally practical to call c this routine with large tables which are not sparse and in c which the 'hybrid' algorithm has little effect. For example, c although it is feasible to compute exact probabilities for the c table c 1 8 5 4 4 2 2 c 5 3 3 4 3 1 0 c 10 1 4 0 0 0 0, c computing exact probabilities for a similar table which has been c enlarged by the addition of an extra row (or column) may not be c feasible. c----------------------------------------------------------------------- subroutine fexact (nrow, ncol, table, ldtabl, expect, percnt, & emin, prt, pre) c SPECIFICATIONS FOR ARGUMENTS integer nrow, ncol, ldtabl double precision expect, percnt, emin, prt, pre, table(ldtabl,*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, i1, i10, i2, i3, i3a, i3b, i3c, i4, i5, i6, i7, & i8, i9, i9a, iiwk, ireal, irwk, iwkmax, iwkpt, & j, k, kk, ldkey, ldstp, mult, nco, nro, & ntot, numb c SPECIFICATIONS FOR INTRINSICS intrinsic max0 integer max0 c SPECIFICATIONS FOR SUBROUTINES external prterr, f2xact c SPECIFICATIONS FOR FUNCTIONS external iwork integer iwork c*********************************************************************** c To increase workspace, increase the c size of of rwrk and set the value of c IWKMAX to the new dimension c c When changing precision, the c following declarations should not be c changed. c*********************************************************************** real rwrk(200000) double precision dwrk(100000) integer iwrk(100000) equivalence (rwrk(1), iwrk(1)), (rwrk(1),dwrk(1)) c Set workspace size iwkmax = 200000 c*********************************************************************** c To increase the length of the table c of paste path lengths relative to the c length of the hash table, increase c MULT c*********************************************************************** mult = 30 c*********************************************************************** c Set IREAL = 4 for DOUBLE PRECISION c Set IREAL = 3 for SINGLE PRECISION c*********************************************************************** ireal = 4 c*********************************************************************** c AMISS is a missing value indicator c which is returned when the c probability is not defined. c*********************************************************************** amiss = -12345.0d0 c iwkpt = 1 c if (nrow .gt. ldtabl) then call prterr (1, 'NROW must be less than or equal to '// & 'LDTABL.') end if ntot = 0 do 20 i=1, nrow do 10 j=1, ncol if (table(i,j) .lt. 0) then call prterr (2, 'All elements of TABLE must '// & 'be positive.') end if ntot = ntot + table(i,j) 10 continue 20 continue if (ntot .eq. 0) then call prterr (3, 'All elements of TABLE are zero. '// & 'PRT and PRE are set to missing values '// & '(NaN, not a number).') prt = amiss pre = amiss go to 9000 end if c nco = max0(nrow,ncol) nro = nrow + ncol - nco k = nrow + ncol + 1 kk = k*max0(nrow,ncol) c i1 = iwork(iwkmax,iwkpt,ntot+1,ireal) i2 = iwork(iwkmax,iwkpt,nco,2) i3 = iwork(iwkmax,iwkpt,nco,2) i3a = iwork(iwkmax,iwkpt,nco,2) i3b = iwork(iwkmax,iwkpt,nro,2) i3c = iwork(iwkmax,iwkpt,nro,2) iiwk = iwork(iwkmax,iwkpt,max0(5*k+2*kk,800+7*max0(nrow,ncol)),2) irwk = iwork(iwkmax,iwkpt,max0(400+max0(nrow,ncol)+1,k),ireal) c Double precision if (ireal .eq. 4) then numb = 18 + 10*mult ldkey = (iwkmax-iwkpt+1)/numb else c Real workspace numb = 12 + 8*mult ldkey = (iwkmax-iwkpt+1)/numb end if c ldstp = mult*ldkey i4 = iwork(iwkmax,iwkpt,2*ldkey,2) i5 = iwork(iwkmax,iwkpt,2*ldkey,2) i6 = iwork(iwkmax,iwkpt,2*ldstp,ireal) i7 = iwork(iwkmax,iwkpt,6*ldstp,2) i8 = iwork(iwkmax,iwkpt,2*ldkey,ireal) i9 = iwork(iwkmax,iwkpt,2*ldkey,ireal) i9a = iwork(iwkmax,iwkpt,2*ldkey,ireal) i10 = iwork(iwkmax,iwkpt,2*ldkey,2) c*********************************************************************** c To convert to double precision, c change RWRK to WWRK in the next CALL c*********************************************************************** c call f2xact (nrow, ncol, table, ldtabl, expect, percnt, emin, & prt, pre, dwrk(i1), iwrk(i2), iwrk(i3), iwrk(i3a), & iwrk(i3b), iwrk(i3c), iwrk(i4), ldkey, iwrk(i5), & dwrk(i6), ldstp, iwrk(i7), dwrk(i8), dwrk(i9), & dwrk(i9a), iwrk(i10), iwrk(iiwk), dwrk(irwk)) c 9000 return end c----------------------------------------------------------------------- c Name: F2XACT c c Purpose: Computes Fisher's exact test for a contingency table, c routine with workspace variables specified. c c Usage: CALL F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT, c EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF, c IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ, c DLP, DSP, TM, KEY2, IWK, RWK) c----------------------------------------------------------------------- subroutine f2xact (nrow, ncol, table, ldtabl, expect, percnt, & emin, prt, pre, fact, ico, iro, kyy, idif, & irn, key, ldkey, ipoin, stp, ldstp, ifrq, & dlp, dsp, tm, key2, iwk, rwk) c SPECIFICATIONS FOR ARGUMENTS integer nrow, ncol, ldtabl, ldkey, ldstp, ico(*), iro(*), & kyy(*), idif(*), irn(*), key(*), ipoin(*), ifrq(*), & key2(*), iwk(*) double precision expect, percnt, emin, prt, pre, table(ldtabl,*), & fact(0:*), stp(*), dlp(*), dsp(*), tm(*), rwk(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, i31, i310, i311, i32, i33, i34, i35, i36, i37, & i38, i39, i41, i42, i43, i44, i45, i46, i47, i48, & iflag, ifreq, ii, ikkey, ikstp, ikstp2, ipn, ipo, & itmp, itop, itp, j, jkey, jstp, jstp2, jstp3, jstp4, & k, k1, kb, kd, kmax, ks, kval, last, n, ncell, nco, & nrb, nro, nro2, ntot, ifault, imax double precision dd, ddf, df, drn, dro, dspt, emn, obs, obs2, & obs3, pastp, pv, tmp, tol logical chisq, ipsh c SPECIFICATIONS FOR INTRINSICS intrinsic dlog, dmax1, dmin1, dexp, max0, min0, mod, nint, dble integer max0, min0, mod, nint double precision dlog, dmax1, dmin1, dexp, dble c SPECIFICATIONS FOR SUBROUTINES external prterr, f3xact, f4xact, f5xact, f6xact, f7xact, isort c SPECIFICATIONS FOR FUNCTIONS external f9xact, gammds double precision f9xact, gammds c*********************************************************************** c IMAX is the largest representable c integer on the machine c*********************************************************************** data imax/2147483647/ c*********************************************************************** c AMISS is a missing value indicator c which is returned when the c probability is not defined. c*********************************************************************** data amiss/-12345.0d0/ c*********************************************************************** c TOL is chosen as the square root of c the smallest relative spacing c*********************************************************************** data tol/3.45254d-07/ c*********************************************************************** c EMX is a large positive value used c in comparing expected values c*********************************************************************** data emx/1.0d30/ c Initialize KEY array do 10 i=1, 2*ldkey key(i) = -9999 key2(i) = -9999 10 continue c Initialize parameters pre = 0.0 itop = 0 if (expect .gt. 0.0d0) then emn = emin else emn = emx end if c Initialize pointers for workspace k = max0(nrow,ncol) c f3xact i31 = 1 i32 = i31 + k i33 = i32 + k i34 = i33 + k i35 = i34 + k i36 = i35 + k i37 = i36 + k i38 = i37 + k i39 = i38 + 400 i310 = 1 i311 = 401 c f4xact k = nrow + ncol + 1 i41 = 1 i42 = i41 + k i43 = i42 + k i44 = i43 + k i45 = i44 + k i46 = i45 + k i47 = i46 + k*max0(nrow,ncol) i48 = 1 c Check table dimensions if (nrow .gt. ldtabl) then call prterr (1, 'NROW must be less than or equal to '// & 'LDTABL.') end if if (ncol .le. 1) then call prterr (4, 'NCOL must be greater than 1.0.') end if c Compute row marginals and total ntot = 0 do 30 i=1, nrow iro(i) = 0 do 20 j=1, ncol if (table(i,j) .lt. -0.0001d0) then call prterr (2, 'All elements of TABLE must be '// & 'positive.') end if iro(i) = iro(i) + nint(table(i,j)) ntot = ntot + nint(table(i,j)) 20 continue 30 continue c if (ntot .eq. 0) then call prterr (3, 'All elements of TABLE are zero. '// & 'PRT and PRE are set to missing values '// & '(NaN, not a number).') prt = amiss pre = amiss go to 9000 end if c Column marginals do 50 i=1, ncol ico(i) = 0 do 40 j=1, nrow ico(i) = ico(i) + nint(table(j,i)) 40 continue 50 continue c sort call isort (nrow, iro) call isort (ncol, ico) c Determine row and column marginals c if (nrow .gt. ncol) then nro = ncol nco = nrow c Interchange row and column marginals do 60 i=1, nrow itmp = iro(i) if (i .le. ncol) iro(i) = ico(i) ico(i) = itmp 60 continue else nro = nrow nco = ncol end if c c Get multiplers for stack kyy(1) = 1 do 70 i=2, nro c Hash table multipliers if (iro(i-1)+1 .le. imax/kyy(i-1)) then kyy(i) = kyy(i-1)*(iro(i-1)+1) j = j/kyy(i-1) else call prterr (5, 'The hash table key cannot be computed'// & ' because the largest key is larger than the'// & ' largest representable integer. The '// & 'algorithm cannot proceed.') end if 70 continue c Maximum product if (iro(nro-1)+1 .le. imax/kyy(nro-1)) then kmax = (iro(nro)+1)*kyy(nro-1) else call prterr (5, 'The hash table key cannot be computed'// & ' because the largest key is larger than the'// & ' largest representable integer. The '// & 'algorithm cannot proceed.') go to 9000 end if c Compute log factorials fact(0) = 0.0d0 fact(1) = 0.0d0 fact(2) = dlog(2.0d0) do 80 i=3, ntot, 2 fact(i) = fact(i-1) + dlog(dble(i)) j = i + 1 if (j .le. ntot) fact(j) = fact(i) + fact(2) + fact(j/2) - & fact(j/2-1) 80 continue c Compute observed path length: OBS obs = tol ntot = 0 do 100 j=1, nco dd = 0.0 do 90 i=1, nro if (nrow .le. ncol) then dd = dd + fact(nint(table(i,j))) ntot = ntot + nint(table(i,j)) else dd = dd + fact(nint(table(j,i))) ntot = ntot + nint(table(j,i)) end if 90 continue obs = obs + fact(ico(j)) - dd 100 continue c Denominator of observed table: DRO dro = f9xact(nro,ntot,iro,fact) prt = dexp(obs-dro) c Initialize pointers k = nco last = ldkey + 1 jkey = ldkey + 1 jstp = ldstp + 1 jstp2 = 3*ldstp + 1 jstp3 = 4*ldstp + 1 jstp4 = 5*ldstp + 1 ikkey = 0 ikstp = 0 ikstp2 = 2*ldstp ipo = 1 ipoin(1) = 1 stp(1) = 0.0 ifrq(1) = 1 ifrq(ikstp2+1) = -1 c 110 kb = nco - k + 1 ks = 0 n = ico(kb) kd = nro + 1 kmax = nro c IDIF is the difference in going to th c daughter do 120 i=1, nro idif(i) = 0 120 continue c Generate the first daughter 130 kd = kd - 1 ntot = min0(n,iro(kd)) idif(kd) = ntot if (idif(kmax) .eq. 0) kmax = kmax - 1 n = n - ntot if (n.gt.0 .and. kd.ne.1) go to 130 if (n .ne. 0) go to 310 c k1 = k - 1 n = ico(kb) ntot = 0 do 140 i=kb + 1, nco ntot = ntot + ico(i) 140 continue c Arc to daughter length=ICO(KB) 150 do 160 i=1, nro irn(i) = iro(i) - idif(i) 160 continue c Sort irn if (k1 .gt. 1) then if (nro .eq. 2) then if (irn(1) .gt. irn(2)) then ii = irn(1) irn(1) = irn(2) irn(2) = ii end if else if (nro .eq. 3) then ii = irn(1) if (ii .gt. irn(3)) then if (ii .gt. irn(2)) then if (irn(2) .gt. irn(3)) then irn(1) = irn(3) irn(3) = ii else irn(1) = irn(2) irn(2) = irn(3) irn(3) = ii end if else irn(1) = irn(3) irn(3) = irn(2) irn(2) = ii end if else if (ii .gt. irn(2)) then irn(1) = irn(2) irn(2) = ii else if (irn(2) .gt. irn(3)) then ii = irn(2) irn(2) = irn(3) irn(3) = ii end if else do 180 j=2, nro i = j - 1 ii = irn(j) 170 if (ii .lt. irn(i)) then irn(i+1) = irn(i) i = i - 1 if (i .gt. 0) go to 170 end if irn(i+1) = ii 180 continue end if c Adjust start for zero do 190 i=1, nro if (irn(i) .ne. 0) go to 200 190 continue 200 nrb = i nro2 = nro - i + 1 else nrb = 1 nro2 = nro end if c Some table values ddf = f9xact(nro,n,idif,fact) drn = f9xact(nro2,ntot,irn(nrb),fact) - dro + ddf c Get hash value if (k1 .gt. 1) then kval = irn(1) + irn(2)*kyy(2) do 210 i=3, nro kval = kval + irn(i)*kyy(i) 210 continue c Get hash table entry i = mod(kval,2*ldkey) + 1 c Search for unused location do 220 itp=i, 2*ldkey ii = key2(itp) if (ii .eq. kval) then go to 240 else if (ii .lt. 0) then key2(itp) = kval dlp(itp) = 1.0d0 dsp(itp) = 1.0d0 go to 240 end if 220 continue c do 230 itp=1, i - 1 ii = key2(itp) if (ii .eq. kval) then go to 240 else if (ii .lt. 0) then key2(itp) = kval dlp(itp) = 1.0 go to 240 end if 230 continue c call prterr (6, 'LDKEY is too small. It is not possible to '// & 'give thevalue of LDKEY required, but you could '// & 'try doubling LDKEY (and possibly LDSTP).') end if c 240 ipsh = .true. c Recover pastp ipn = ipoin(ipo+ikkey) pastp = stp(ipn+ikstp) ifreq = ifrq(ipn+ikstp) c Compute shortest and longest path if (k1 .gt. 1) then obs2 = obs - fact(ico(kb+1)) - fact(ico(kb+2)) - ddf do 250 i=3, k1 obs2 = obs2 - fact(ico(kb+i)) 250 continue c if (dlp(itp) .gt. 0.0d0) then dspt = obs - obs2 - ddf c Compute longest path dlp(itp) = 0.0d0 call f3xact (nro2, irn(nrb), k1, ico(kb+1), dlp(itp), & ntot, fact, iwk(i31), iwk(i32), iwk(i33), & iwk(i34), iwk(i35), iwk(i36), iwk(i37), & iwk(i38), iwk(i39), rwk(i310), rwk(i311), tol) dlp(itp) = dmin1(0.0d0,dlp(itp)) c Compute shortest path dsp(itp) = dspt call f4xact (nro2, irn(nrb), k1, ico(kb+1), dsp(itp), & fact, iwk(i47), iwk(i41), iwk(i42), iwk(i43), & iwk(i44), iwk(i45), iwk(i46), rwk(i48), tol) dsp(itp) = dmin1(0.0d0,dsp(itp)-dspt) c Use chi-squared approximation? if (dble(irn(nrb)*ico(kb+1))/dble(ntot) .gt. emn) then ncell = 0.0 do 270 i=1, nro2 do 260 j=1, k1 if (irn(nrb+i-1)*ico(kb+j) .ge. ntot*expect) then ncell = ncell + 1 end if 260 continue 270 continue if (ncell*100 .ge. k1*nro2*percnt) then tmp = 0.0 do 280 i=1, nro2 tmp = tmp + fact(irn(nrb+i-1)) - & fact(irn(nrb+i-1)-1) 280 continue tmp = tmp*(k1-1) do 290 j=1, k1 tmp = tmp + (nro2-1)*(fact(ico(kb+j))-fact(ico(kb+ & j)-1)) 290 continue df = (nro2-1)*(k1-1) tmp = tmp + df*1.83787706640934548356065947281d0 tmp = tmp - (nro2*k1-1)*(fact(ntot)-fact(ntot-1)) tm(itp) = -2.0d0*(obs-dro) - tmp else c tm(itp) set to a flag value tm(itp) = -9876.0d0 end if else tm(itp) = -9876.0d0 end if end if obs3 = obs2 - dlp(itp) obs2 = obs2 - dsp(itp) if (tm(itp) .eq. -9876.0d0) then chisq = .false. else chisq = .true. tmp = tm(itp) end if else obs2 = obs - drn - dro obs3 = obs2 end if c Process node with new PASTP 300 if (pastp .le. obs3) then c Update pre pre = pre + dble(ifreq)*dexp(pastp+drn) c else if (pastp .lt. obs2) then if (chisq) then df = (nro2-1)*(k1-1) pv = 1.0 - gammds(dmax1(0.0d0,tmp+2.0d0*(pastp+drn))/ & 2.0d0,df/2.0d0,ifault) pre = pre + dble(ifreq)*dexp(pastp+drn)*pv else c Put daughter on queue call f5xact (pastp+ddf, tol, kval, key(jkey), ldkey, & ipoin(jkey), stp(jstp), ldstp, ifrq(jstp), & ifrq(jstp2), ifrq(jstp3), ifrq(jstp4), ifreq, & itop, ipsh) ipsh = .false. end if end if c Get next PASTP on chain ipn = ifrq(ipn+ikstp2) if (ipn .gt. 0) then pastp = stp(ipn+ikstp) ifreq = ifrq(ipn+ikstp) go to 300 end if c Generate a new daughter node call f7xact (kmax, iro, idif, kd, ks, iflag) if (iflag .ne. 1) go to 150 c Go get a new mother from stage K 310 iflag = 1 call f6xact (nro, iro, iflag, kyy, key(ikkey+1), ldkey, last, & ipo) c Update pointers if (iflag .eq. 3) then k = k - 1 itop = 0 ikkey = jkey - 1 ikstp = jstp - 1 ikstp2 = jstp2 - 1 jkey = ldkey - jkey + 2 jstp = ldstp - jstp + 2 jstp2 = 2*ldstp + jstp do 320 i=1, 2*ldkey key2(i) = -9999 320 continue if (k .ge. 2) go to 310 else go to 110 end if c 9000 return end c----------------------------------------------------------------------- c Name: F3XACT c c Purpose: Computes the shortest path length for a given table. c c Usage: CALL F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, c IRO, IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, c TOL) c c Arguments: c NROW - The number of rows in the table. (Input) c IROW - Vector of length NROW containing the row sums for the c table. (Input) c NCOL - The number of columns in the table. (Input) c ICOL - Vector of length K containing the column sums for the c table. (Input) c DLP - The longest path for the table. (Output) c MM - The total count in the table. (Output) c FACT - Vector containing the logarithms of factorials. (Input) c ICO - Work vector of length MAX(NROW,NCOL). c IRO - Work vector of length MAX(NROW,NCOL). c IT - Work vector of length MAX(NROW,NCOL). c LB - Work vector of length MAX(NROW,NCOL). c NR - Work vector of length MAX(NROW,NCOL). c NT - Work vector of length MAX(NROW,NCOL). c NU - Work vector of length MAX(NROW,NCOL). c ITC - Work vector of length 400. c IST - Work vector of length 400. c STV - Work vector of length 400. c ALEN - Work vector of length MAX(NROW,NCOL). c TOL - Tolerance. (Input) c----------------------------------------------------------------------- subroutine f3xact (nrow, irow, ncol, icol, dlp, mm, fact, ico, & iro, it, lb, nr, nt, nu, itc, ist, stv, alen, & tol) c SPECIFICATIONS FOR ARGUMENTS integer nrow, ncol, mm, irow(*), icol(*), ico(*), iro(*), & it(*), lb(*), nr(*), nt(*), nu(*), itc(*), ist(*) double precision dlp, tol, fact(0:*), stv(*), alen(0:*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, ic1, ic2, ii, ipn, irl, itp, k, key, ks, kyy, lev, & n11, n12, nc1, nc1s, nco, nct, nn, nn1, nr1, nro, nrt double precision v, val, vmn logical xmin c SPECIFICATIONS FOR SAVE VARIABLES integer ldst, nitc, nst save ldst, nitc, nst c SPECIFICATIONS FOR INTRINSICS intrinsic dmin1, int, mod, dble integer int, mod double precision dmin1, dble c SPECIFICATIONS FOR SUBROUTINES external prterr, f10act, isort c data ldst/200/, nst/0/, nitc/0/ c do 10 i=0, ncol alen(i) = 0.0 10 continue do 20 i=1, 400 ist(i) = -1 20 continue c nrow is 1 if (nrow .le. 1) then if (nrow .gt. 0) then dlp = dlp - fact(icol(1)) do 30 i=2, ncol dlp = dlp - fact(icol(i)) 30 continue end if go to 9000 end if c ncol is 1 if (ncol .le. 1) then if (ncol .gt. 0) then dlp = dlp - fact(irow(1)) - fact(irow(2)) do 40 i=3, nrow dlp = dlp - fact(irow(i)) 40 continue end if go to 9000 end if c 2 by 2 table if (nrow*ncol .eq. 4) then n11 = (irow(1)+1)*(icol(1)+1)/(mm+2) n12 = irow(1) - n11 dlp = dlp - fact(n11) - fact(n12) - fact(icol(1)-n11) - & fact(icol(2)-n12) go to 9000 end if c Test for optimal table val = 0.0 xmin = .false. if (irow(nrow) .le. irow(1)+ncol) then call f10act (nrow, irow, ncol, icol, val, xmin, fact, lb, nu, & nr) end if if (.not.xmin) then if (icol(ncol) .le. icol(1)+nrow) then call f10act (ncol, icol, nrow, irow, val, xmin, fact, lb, & nu, nr) end if end if c if (xmin) then dlp = dlp - val go to 9000 end if c Setup for dynamic programming nn = mm c Minimize ncol if (nrow .ge. ncol) then nro = nrow nco = ncol c do 50 i=1, nrow iro(i) = irow(i) 50 continue c ico(1) = icol(1) nt(1) = nn - ico(1) do 60 i=2, ncol ico(i) = icol(i) nt(i) = nt(i-1) - ico(i) 60 continue else nro = ncol nco = nrow c ico(1) = irow(1) nt(1) = nn - ico(1) do 70 i=2, nrow ico(i) = irow(i) nt(i) = nt(i-1) - ico(i) 70 continue c do 80 i=1, ncol iro(i) = icol(i) 80 continue end if c Initialize pointers vmn = 1.0d10 nc1s = nco - 1 irl = 1 ks = 0 k = ldst kyy = ico(nco) + 1 go to 100 c Test for optimality 90 xmin = .false. if (iro(nro) .le. iro(irl)+nco) then call f10act (nro, iro(irl), nco, ico, val, xmin, fact, lb, & nu, nr) end if if (.not.xmin) then if (ico(nco) .le. ico(1)+nro) then call f10act (nco, ico, nro, iro(irl), val, xmin, fact, lb, & nu, nr) end if end if c if (xmin) then if (val .lt. vmn) vmn = val go to 200 end if c Setup to generate new node 100 lev = 1 nr1 = nro - 1 nrt = iro(irl) nct = ico(1) lb(1) = int(dble((nrt+1)*(nct+1))/dble(nn+nr1*nc1s+1)-tol) - 1 nu(1) = int(dble((nrt+nc1s)*(nct+nr1))/dble(nn+nr1+nc1s)) - & lb(1) + 1 nr(1) = nrt - lb(1) c Generate a node 110 nu(lev) = nu(lev) - 1 if (nu(lev) .eq. 0) then if (lev .eq. 1) go to 200 lev = lev - 1 go to 110 end if lb(lev) = lb(lev) + 1 nr(lev) = nr(lev) - 1 120 alen(lev) = alen(lev-1) + fact(lb(lev)) if (lev .lt. nc1s) then nn1 = nt(lev) nrt = nr(lev) lev = lev + 1 nc1 = nco - lev nct = ico(lev) lb(lev) = dble((nrt+1)*(nct+1))/dble(nn1+nr1*nc1+1) - tol nu(lev) = dble((nrt+nc1)*(nct+nr1))/dble(nn1+nr1+nc1) - & lb(lev) + 1 nr(lev) = nrt - lb(lev) go to 120 end if alen(nco) = alen(lev) + fact(nr(lev)) lb(nco) = nr(lev) c v = val + alen(nco) if (nro .eq. 2) then c Only 1 row left v = v + fact(ico(1)-lb(1)) + fact(ico(2)-lb(2)) do 130 i=3, nco v = v + fact(ico(i)-lb(i)) 130 continue if (v .lt. vmn) vmn = v else if (nro.eq.3 .and. nco.eq.2) then c 3 rows and 2 columns nn1 = nn - iro(irl) + 2 ic1 = ico(1) - lb(1) ic2 = ico(2) - lb(2) n11 = (iro(irl+1)+1)*(ic1+1)/nn1 n12 = iro(irl+1) - n11 v = v + fact(n11) + fact(n12) + fact(ic1-n11) + & fact(ic2-n12) if (v .lt. vmn) vmn = v else c Column marginals are new node do 140 i=1, nco it(i) = ico(i) - lb(i) 140 continue c Sort column marginals if (nco .eq. 2) then if (it(1) .gt. it(2)) then ii = it(1) it(1) = it(2) it(2) = ii end if else if (nco .eq. 3) then ii = it(1) if (ii .gt. it(3)) then if (ii .gt. it(2)) then if (it(2) .gt. it(3)) then it(1) = it(3) it(3) = ii else it(1) = it(2) it(2) = it(3) it(3) = ii end if else it(1) = it(3) it(3) = it(2) it(2) = ii end if else if (ii .gt. it(2)) then it(1) = it(2) it(2) = ii else if (it(2) .gt. it(3)) then ii = it(2) it(2) = it(3) it(3) = ii end if else call isort (nco, it) end if c Compute hash value key = it(1)*kyy + it(2) do 150 i=3, nco key = it(i) + key*kyy 150 continue c Table index ipn = mod(key,ldst) + 1 c Find empty position ii = ks + ipn do 160 itp=ipn, ldst if (ist(ii) .lt. 0) then go to 180 else if (ist(ii) .eq. key) then go to 190 end if ii = ii + 1 160 continue c ii = ks + 1 do 170 itp=1, ipn - 1 if (ist(ii) .lt. 0) then go to 180 else if (ist(ii) .eq. key) then go to 190 end if ii = ii + 1 170 continue c call prterr (30, 'Stack length exceeded in f3xact.'// & ' This problem should not occur.') c Push onto stack 180 ist(ii) = key stv(ii) = v nst = nst + 1 ii = nst + ks itc(ii) = itp go to 110 c Marginals already on stack 190 stv(ii) = dmin1(v,stv(ii)) end if go to 110 c Pop item from stack 200 if (nitc .gt. 0) then c Stack index itp = itc(nitc+k) + k nitc = nitc - 1 val = stv(itp) key = ist(itp) ist(itp) = -1 c Compute marginals do 210 i=nco, 2, -1 ico(i) = mod(key,kyy) key = key/kyy 210 continue ico(1) = key c Set up nt array nt(1) = nn - ico(1) do 220 i=2, nco nt(i) = nt(i-1) - ico(i) 220 continue go to 90 c else if (nro.gt.2 .and. nst.gt.0) then c Go to next level nitc = nst nst = 0 k = ks ks = ldst - ks nn = nn - iro(irl) irl = irl + 1 nro = nro - 1 go to 200 end if c dlp = dlp - vmn 9000 return end c----------------------------------------------------------------------- c Name: F4XACT c c Purpose: Computes the longest path length for a given table. c c Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK, c NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK, c TOL) c c Arguments: c NROW - The number of rows in the table. (Input) c IROW - Vector of length NROW containing the row sums for the c table. (Input) c NCOL - The number of columns in the table. (Input) c ICOL - Vector of length K containing the column sums for the c table. (Input) c DSP - The shortest path for the table. (Output) c FACT - Vector containing the logarithms of factorials. (Input) c ICSTK - NCOL by NROW+NCOL+1 work array. c NCSTK - Work vector of length NROW+NCOL+1. c LSTK - Work vector of length NROW+NCOL+1. c MSTK - Work vector of length NROW+NCOL+1. c NSTK - Work vector of length NROW+NCOL+1. c NRSTK - Work vector of length NROW+NCOL+1. c IRSTK - NROW by MAX(NROW,NCOL) work array. c YSTK - Work vector of length NROW+NCOL+1. c TOL - Tolerance. (Input) c----------------------------------------------------------------------- subroutine f4xact (nrow, irow, ncol, icol, dsp, fact, icstk, & ncstk, lstk, mstk, nstk, nrstk, irstk, ystk, & tol) c SPECIFICATIONS FOR ARGUMENTS integer nrow, ncol, irow(*), icol(*), icstk(ncol,*), & ncstk(*), lstk(*), mstk(*), nstk(*), nrstk(*), & irstk(nrow,*) double precision dsp, tol, fact(0:*), ystk(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, ic1, ict, ir1, irt, istk, j, k, l, m, mn, n, nco, & nro double precision amx, y c SPECIFICATIONS FOR SUBROUTINES external f11act, f8xact c Take care of the easy cases firstkt if (nrow .eq. 1) then do 10 i=1, ncol dsp = dsp - fact(icol(i)) 10 continue go to 9000 end if c if (ncol .eq. 1) then do 20 i=1, nrow dsp = dsp - fact(irow(i)) 20 continue go to 9000 end if c if (nrow*ncol .eq. 4) then if (irow(2) .le. icol(2)) then dsp = dsp - fact(irow(2)) - fact(icol(1)) - & fact(icol(2)-irow(2)) else dsp = dsp - fact(icol(2)) - fact(irow(1)) - & fact(irow(2)-icol(2)) end if go to 9000 end if c initialization before loop do 30 i=1, nrow irstk(i,1) = irow(nrow-i+1) 30 continue c do 40 j=1, ncol icstk(j,1) = icol(ncol-j+1) 40 continue c nro = nrow nco = ncol nrstk(1) = nro ncstk(1) = nco ystk(1) = 0.0 y = 0.0 istk = 1 l = 1 amx = 0.0 c 50 ir1 = irstk(1,istk) ic1 = icstk(1,istk) if (ir1 .gt. ic1) then if (nro .ge. nco) then m = nco - 1 n = 2 else m = nro n = 1 end if else if (ir1 .lt. ic1) then if (nro .le. nco) then m = nro - 1 n = 1 else m = nco n = 2 end if else if (nro .le. nco) then m = nro - 1 n = 1 else m = nco - 1 n = 2 end if end if c 60 if (n .eq. 1) then i = l j = 1 else i = 1 j = l end if c irt = irstk(i,istk) ict = icstk(j,istk) mn = irt if (mn .gt. ict) mn = ict y = y + fact(mn) if (irt .eq. ict) then nro = nro - 1 nco = nco - 1 call f11act (irstk(1,istk), i, nro, irstk(1,istk+1)) call f11act (icstk(1,istk), j, nco, icstk(1,istk+1)) else if (irt .gt. ict) then nco = nco - 1 call f11act (icstk(1,istk), j, nco, icstk(1,istk+1)) call f8xact (irstk(1,istk), irt-ict, i, nro, irstk(1,istk+1)) else nro = nro - 1 call f11act (irstk(1,istk), i, nro, irstk(1,istk+1)) call f8xact (icstk(1,istk), ict-irt, j, nco, icstk(1,istk+1)) end if c if (nro .eq. 1) then do 70 k=1, nco y = y + fact(icstk(k,istk+1)) 70 continue go to 90 end if c if (nco .eq. 1) then do 80 k=1, nro y = y + fact(irstk(k,istk+1)) 80 continue go to 90 end if c lstk(istk) = l mstk(istk) = m nstk(istk) = n istk = istk + 1 nrstk(istk) = nro ncstk(istk) = nco ystk(istk) = y l = 1 go to 50 c 90 if (y .gt. amx) then amx = y if (dsp-amx .le. tol) then dsp = 0.0 go to 9000 end if end if c 100 istk = istk - 1 if (istk .eq. 0) then dsp = dsp - amx if (dsp-amx .le. tol) dsp = 0.0 go to 9000 end if l = lstk(istk) + 1 c 110 if (l .gt. mstk(istk)) go to 100 n = nstk(istk) nro = nrstk(istk) nco = ncstk(istk) y = ystk(istk) if (n .eq. 1) then if (irstk(l,istk) .lt. irstk(l-1,istk)) go to 60 else if (n .eq. 2) then if (icstk(l,istk) .lt. icstk(l-1,istk)) go to 60 end if c l = l + 1 go to 110 9000 return end c----------------------------------------------------------------------- c Name: F5XACT c c Purpose: Put node on stack in network algorithm. c c Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP, c LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP, c IPSH) c c Arguments: c PASTP - The past path length. (Input) c TOL - Tolerance for equivalence of past path lengths. (Input) c KVAL - Key value. (Input) c KEY - Vector of length LDKEY containing the key values. c (Input/output) c LDKEY - Length of vector KEY. (Input) c IPOIN - Vector of length LDKEY pointing to the linked list c of past path lengths. (Input/output) c STP - Vector of length LSDTP containing the linked lists c of past path lengths. (Input/output) c LDSTP - Length of vector STP. (Input) c IFRQ - Vector of length LDSTP containing the past path c frequencies. (Input/output) c NPOIN - Vector of length LDSTP containing the pointers to c the next past path length. (Input/output) c NR - Vector of length LDSTP containing the right object c pointers in the tree of past path lengths. c (Input/output) c NL - Vector of length LDSTP containing the left object c pointers in the tree of past path lengths. c (Input/output) c IFREQ - Frequency of the current path length. (Input) c ITOP - Pointer to the top of STP. (Input) c IPSH - Option parameter. (Input) c If IPSH is true, the past path length is found in the c table KEY. Otherwise the location of the past path c length is assumed known and to have been found in c a previous call. c----------------------------------------------------------------------- subroutine f5xact (pastp, tol, kval, key, ldkey, ipoin, stp, & ldstp, ifrq, npoin, nr, nl, ifreq, itop, ipsh) c SPECIFICATIONS FOR ARGUMENTS integer kval, ldkey, ldstp, ifreq, itop, key(*), ipoin(*), & ifrq(*), npoin(*), nr(*), nl(*) double precision pastp, tol, stp(*) logical ipsh c SPECIFICATIONS FOR LOCAL VARIABLES integer ipn, ird, itmp double precision test1, test2 c SPECIFICATIONS FOR SAVE VARIABLES integer itp save itp c SPECIFICATIONS FOR INTRINSICS intrinsic mod integer mod c SPECIFICATIONS FOR SUBROUTINES external prterr c if (ipsh) then c Convert KVAL to integer in range c 1, ..., LDKEY. ird = mod(kval,ldkey) + 1 c Search for an unused location do 10 itp=ird, ldkey if (key(itp) .eq. kval) go to 40 if (key(itp) .lt. 0) go to 30 10 continue c do 20 itp=1, ird - 1 if (key(itp) .eq. kval) go to 40 if (key(itp) .lt. 0) go to 30 20 continue c Return if KEY array is full call prterr(6, 'LDKEY is too small for this problem. It is '// & 'not possible to estimate the value of LDKEY '// & 'required, but twice the current value may be '// & 'sufficient.') c Update KEY 30 key(itp) = kval itop = itop + 1 ipoin(itp) = itop c Return if STP array full if (itop .gt. ldstp) then call prterr(7, 'LDSTP is too small for this problem. It '// & 'is not possible to estimate the value of '// & 'LDSTP required, but twice the current value '// & 'may be sufficient.') end if c Update STP, etc. npoin(itop) = -1 nr(itop) = -1 nl(itop) = -1 stp(itop) = pastp ifrq(itop) = ifreq go to 9000 end if c Find location, if any, of pastp 40 ipn = ipoin(itp) test1 = pastp - tol test2 = pastp + tol c 50 if (stp(ipn) .lt. test1) then ipn = nl(ipn) if (ipn .gt. 0) go to 50 else if (stp(ipn) .gt. test2) then ipn = nr(ipn) if (ipn .gt. 0) go to 50 else ifrq(ipn) = ifrq(ipn) + ifreq go to 9000 end if c Return if STP array full itop = itop + 1 if (itop .gt. ldstp) then call prterr(7, 'LDSTP is too small for this problem. It is '// & 'not possible to estimate the value of LDSTP '// & 'rerquired, but twice the current value may be '// & 'sufficient.') go to 9000 end if c Find location to add value ipn = ipoin(itp) itmp = ipn 60 if (stp(ipn) .lt. test1) then itmp = ipn ipn = nl(ipn) if (ipn .gt. 0) then go to 60 else nl(itmp) = itop end if else if (stp(ipn) .gt. test2) then itmp = ipn ipn = nr(ipn) if (ipn .gt. 0) then go to 60 else nr(itmp) = itop end if end if c Update STP, etc. npoin(itop) = npoin(itmp) npoin(itmp) = itop stp(itop) = pastp ifrq(itop) = ifreq nl(itop) = -1 nr(itop) = -1 c 9000 return end c----------------------------------------------------------------------- c Name: F6XACT c c Purpose: Pop a node off the stack. c c Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, c IPN) c c Arguments: c NROW - The number of rows in the table. (Input) c IROW - Vector of length nrow containing the row sums on output. c (Output) c IFLAG - Set to 3 if there are no additional nodes to process. c (Output) c KYY - Constant mutlipliers used in forming the hash table key. c (Input) c KEY - Vector of length LDKEY containing the hash table keys. c (Input/output) c LDKEY - Length of vector KEY. (Input) c LAST - Index of the last key popped off the stack. c (Input/output) c IPN - Pointer to the linked list of past path lengths. c (Output) c----------------------------------------------------------------------- subroutine f6xact (nrow, irow, iflag, kyy, key, ldkey, last, ipn) c SPECIFICATIONS FOR ARGUMENTS integer nrow, iflag, ldkey, last, ipn, irow(*), kyy(*), key(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer j, kval c SPECIFICATIONS FOR SAVE VARIABLES c 10 last = last + 1 if (last .le. ldkey) then if (key(last) .lt. 0) go to 10 c Get KVAL from the stack kval = key(last) key(last) = -9999 do 20 j=nrow, 2, -1 irow(j) = kval/kyy(j) kval = kval - irow(j)*kyy(j) 20 continue irow(1) = kval ipn = last else last = 0 iflag = 3 end if return end c----------------------------------------------------------------------- c Name: F7XACT c c Purpose: Generate the new nodes for given marinal totals. c c Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG) c c Arguments: c NROW - The number of rows in the table. (Input) c IMAX - The row marginal totals. (Input) c IDIF - The column counts for the new column. (Input/output) c K - Indicator for the row to decrement. (Input/output) c KS - Indicator for the row to increment. (Input/output) c IFLAG - Status indicator. (Output) c If IFLAG is zero, a new table was generated. For c IFLAG = 1, no additional tables could be generated. c----------------------------------------------------------------------- subroutine f7xact (nrow, imax, idif, k, ks, iflag) c SPECIFICATIONS FOR ARGUMENTS integer nrow, k, ks, iflag, imax(*), idif(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, k1, m, mm c SPECIFICATIONS FOR INTRINSICS intrinsic min0 integer min0 c iflag = 0 c Find node which can be c incremented, ks if (ks .eq. 0) then 10 ks = ks + 1 if (idif(ks) .eq. imax(ks)) go to 10 end if c Find node to decrement (>ks) 20 if (idif(k).gt.0 .and. k.gt.ks) then idif(k) = idif(k) - 1 30 k = k - 1 if (imax(k) .eq. 0) go to 30 m = k c Find node to increment (>=ks) 40 if (idif(m) .ge. imax(m)) then m = m - 1 go to 40 end if idif(m) = idif(m) + 1 c Change ks if (m .eq. ks) then if (idif(m) .eq. imax(m)) ks = k end if else c Check for finish 50 do 60 k1=k + 1, nrow if (idif(k1) .gt. 0) go to 70 60 continue iflag = 1 go to 9000 c Reallocate counts 70 mm = 1 do 80 i=1, k mm = mm + idif(i) idif(i) = 0 80 continue k = k1 90 k = k - 1 m = min0(mm,imax(k)) idif(k) = m mm = mm - m if (mm.gt.0 .and. k.ne.1) go to 90 c Check that all counts c reallocated if (mm .gt. 0) then if (k1 .ne. nrow) then k = k1 go to 50 end if iflag = 1 go to 9000 end if c Get ks idif(k1) = idif(k1) - 1 ks = 0 100 ks = ks + 1 if (ks .gt. k) go to 9000 if (idif(ks) .ge. imax(ks)) go to 100 end if c 9000 return end c----------------------------------------------------------------------- c Name: F8XACT c c Purpose: Routine for reducing a vector when there is a zero c element. c c Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW) c c Arguments: c IROW - Vector containing the row counts. (Input) c IS - Indicator. (Input) c I1 - Indicator. (Input) c IZERO - Position of the zero. (Input) c NEW - Vector of new row counts. (Output) c----------------------------------------------------------------------- subroutine f8xact (irow, is, i1, izero, new) c SPECIFICATIONS FOR ARGUMENTS integer is, i1, izero, irow(*), new(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i c do 10 i=1, i1 - 1 new(i) = irow(i) 10 continue c do 20 i=i1, izero - 1 if (is .ge. irow(i+1)) go to 30 new(i) = irow(i+1) 20 continue c i = izero 30 new(i) = is 40 i = i + 1 if (i .gt. izero) return new(i) = irow(i) go to 40 end c----------------------------------------------------------------------- c Name: F9XACT c c Purpose: Computes the log of a multinomial coefficient. c c Usage: F9XACT(N, MM, IR, FACT) c c Arguments: c N - Length of IR. (Input) c MM - Number for factorial in numerator. (Input) c IR - Vector of length N containing the numebers for the c denominator of the factorial. (Input) c FACT - Table of log factorials. (Input) c F9XACT - The log of the multinomal coefficient. (Output) c----------------------------------------------------------------------- double precision function f9xact (n, mm, ir, fact) c SPECIFICATIONS FOR ARGUMENTS integer n, mm, ir(*) double precision fact(0:*) c SPECIFICATIONS FOR LOCAL VARIABLES integer k c f9xact = fact(mm) do 10 k=1, n f9xact = f9xact - fact(ir(k)) 10 continue c return end c----------------------------------------------------------------------- c Name: F10ACT c c Purpose: Computes the shortest path length for special tables. c c Usage: CALL F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, c NE, M) c c Arguments: c NROW - The number of rows in the table. (Input) c IROW - Vector of length NROW containing the row totals. (Input) c NCOL - The number of columns in the table. (Input) c ICO - Vector of length NCOL containing the column totals. c (Input) c VAL - The shortest path. (Output) c XMIN - Set to true if shortest path obtained. (Output) c FACT - Vector containing the logarithms of factorials. c (Input) c ND - Workspace vector of length NROW. c NE - Workspace vector of length NCOL. c M - Workspace vector of length NCOL. c c Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis c----------------------------------------------------------------------- subroutine f10act (nrow, irow, ncol, icol, val, xmin, fact, nd, & ne, m) c SPECIFICATIONS FOR ARGUMENTS integer nrow, ncol, irow(*), icol(*), nd(*), ne(*), m(*) double precision val, fact(0:*) logical xmin c SPECIFICATIONS FOR LOCAL VARIABLES integer i, is, ix, nrw1 c do 10 i=1, nrow - 1 nd(i) = 0 10 continue c is = icol(1)/nrow ne(1) = is ix = icol(1) - nrow*is m(1) = ix if (ix .ne. 0) nd(ix) = nd(ix) + 1 c do 20 i=2, ncol ix = icol(i)/nrow ne(i) = ix is = is + ix ix = icol(i) - nrow*ix m(i) = ix if (ix .ne. 0) nd(ix) = nd(ix) + 1 20 continue c do 30 i=nrow - 2, 1, -1 nd(i) = nd(i) + nd(i+1) 30 continue c ix = 0 nrw1 = nrow + 1 do 40 i=nrow, 2, -1 ix = ix + is + nd(nrw1-i) - irow(i) if (ix .lt. 0) return 40 continue c do 50 i=1, ncol ix = ne(i) is = m(i) val = val + is*fact(ix+1) + (nrow-is)*fact(ix) 50 continue xmin = .true. c return end c----------------------------------------------------------------------- c Name: F11ACT c c Purpose: Routine for revising row totals. c c Usage: CALL F11ACT (IROW, I1, I2, NEW) c c Arguments: c IROW - Vector containing the row totals. (Input) c I1 - Indicator. (Input) c I2 - Indicator. (Input) c NEW - Vector containing the row totals. (Input) c----------------------------------------------------------------------- subroutine f11act (irow, i1, i2, new) c SPECIFICATIONS FOR ARGUMENTS integer i1, i2, irow(*), new(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i c do 10 i=1, i1 - 1 new(i) = irow(i) 10 continue c do 20 i=i1, i2 new(i) = irow(i+1) 20 continue c return end c----------------------------------------------------------------------- c Name: ERPRT c c Purpose: Print an error message and stop. c c Usage: CALL ERPRT (ICODE, MES) c c Arguments: c ICODE - Integer code for the error message. (Input) c MES - Character string containing the error message. (Input) c----------------------------------------------------------------------- subroutine prterr (icode, mes) c SPECIFICATIONS FOR ARGUMENTS integer icode character mes*(*) c write (*,*) 'FEXACT ERROR: ', icode, ' ', mes stop end c----------------------------------------------------------------------- c Name: IWORK c c Purpose: Routine for allocating workspace. c c Usage: IWORK (IWKMAX, IWKPT, NUMBER, ITYPE) c c Arguments: c IWKMAX - Maximum length of workspace. (Input) c IWKPT - Amount of workspace currently allocated. (Input/output) c NUMBER - Number of elements of workspace desired. (Input) c ITYPE - Worspace type. (Input) c ITYPE TYPE c 2 Integer c 3 Real c 4 Double Precision c IWORK - Index in RWRK, DWRK, or IWRK of the beginning of the c first element in the workspace array. (Output) c----------------------------------------------------------------------- integer function iwork (iwkmax, iwkpt, number, itype) c SPECIFICATIONS FOR ARGUMENTS integer iwkmax, iwkpt, number, itype c SPECIFICATIONS FOR INTRINSICS intrinsic mod integer mod c SPECIFICATIONS FOR SUBROUTINES external prterr c iwork = iwkpt if (itype.eq.2 .or. itype.eq.3) then iwkpt = iwkpt + number else if (mod(iwork,2) .ne. 0) iwork = iwork + 1 iwkpt = iwkpt + 2*number iwork = iwork/2 end if if (iwkpt .gt. iwkmax+1) then call prterr (40, 'Out of workspace.') end if return end c----------------------------------------------------------------------- c Name: ISORT c c Purpose: Shell sort for an integer vector. c c Usage: CALL ISORT (N, IX) c c Arguments: c N - Lenth of vector IX. (Input) c IX - Vector to be sorted. (Input/output) c----------------------------------------------------------------------- subroutine isort (n, ix) c SPECIFICATIONS FOR ARGUMENTS integer n, ix(*) c SPECIFICATIONS FOR LOCAL VARIABLES integer i, ikey, il(10), it, iu(10), j, kl, ku, m c SPECIFICATIONS FOR SUBROUTINES external prterr c Sort IX m = 1 i = 1 j = n 10 if (i .ge. j) go to 40 kl = i ku = j ikey = i j = j + 1 c Find element in first half 20 i = i + 1 if (i .lt. j) then if (ix(ikey) .gt. ix(i)) go to 20 end if c Find element in second half 30 j = j - 1 if (ix(j) .gt. ix(ikey)) go to 30 c Interchange if (i .lt. j) then it = ix(i) ix(i) = ix(j) ix(j) = it go to 20 end if it = ix(ikey) ix(ikey) = ix(j) ix(j) = it c Save upper and lower subscripts of c the array yet to be sorted if (m .lt. 11) then if (j-kl .lt. ku-j) then il(m) = j + 1 iu(m) = ku i = kl j = j - 1 else il(m) = kl iu(m) = j - 1 i = j + 1 j = ku end if m = m + 1 go to 10 else call prterr (20, 'This should never occur.') end if c Use another segment 40 m = m - 1 if (m .eq. 0) go to 9000 i = il(m) j = iu(m) go to 10 c 9000 return end c----------------------------------------------------------------------- c Name: GAMMDS c c Purpose: Cumulative distribution for the gamma distribution. c c Usage: PGAMMA (Q, ALPHA,IFAULT) c c Arguments: c Q - Value at which the distribution is desired. (Input) c ALPHA - Parameter in the gamma distribution. (Input) c IFAULT - Error indicator. (Output) c IFAULT DEFINITION c 0 No error c 1 An argument is misspecified. c 2 A numerical error has occurred. c PGAMMA - The cdf for the gamma distribution with parameter alpha c evaluated at Q. (Output) c----------------------------------------------------------------------- c c Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113 c c Computes the incomplete gamma integral for positive c parameters Y, P using and infinite series. c double precision function gammds (y, p, ifault) c SPECIFICATIONS FOR ARGUMENTS integer ifault double precision y, p c SPECIFICATIONS FOR LOCAL VARIABLES integer ifail double precision a, c, f c SPECIFICATIONS FOR SAVE VARIABLES double precision e, one, zero save e, one, zero c SPECIFICATIONS FOR INTRINSICS intrinsic dlog, dexp double precision dlog, dexp c SPECIFICATIONS FOR FUNCTIONS external alogam double precision alogam double precision zexp, zlog c data e, zero, one/1.0d-6, 0.0d0, 1.0d0/ c zexp(a) = dexp(a) zlog(a) = dlog(a) c c Checks for the admissibility of arguments and value of F c ifault = 1 gammds = zero if (y.le.zero .or. p.le.zero) return ifault = 2 c c ALOGAM is natural log of gamma function c no need to test ifail as an error is impossible c f = zexp(p*zlog(y)-alogam(p+one,ifail)-y) if (f .eq. zero) return ifault = 0 c c Series begins c c = one gammds = one a = p 10 a = a + one c = c*y/a gammds = gammds + c if (c/gammds .gt. e) go to 10 gammds = gammds*f return end c----------------------------------------------------------------------- c Name: ALOGAM c c Purpose: Value of the log-gamma function. c c Usage: ALOGAM (X, IFAULT) c c Arguments: c X - Value at which the log-gamma function is to be evaluated. c (Input) c IFAULT - Error indicator. (Output) c IFAULT DEFINITION c 0 No error c 1 X .LT. 0 c ALGAMA - The value of the log-gamma function at XX. (Output) c----------------------------------------------------------------------- c c Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684 c c Evaluates natural logarithm of gamma(x) c for X greater than zero. c double precision function alogam (x, ifault) c SPECIFICATIONS FOR ARGUMENTS integer ifault double precision x c SPECIFICATIONS FOR LOCAL VARIABLES double precision f, y, z c SPECIFICATIONS FOR SAVE VARIABLES double precision a1, a2, a3, a4, a5, half, one, seven, zero save a1, a2, a3, a4, a5, half, one, seven, zero c SPECIFICATIONS FOR INTRINSICS intrinsic dlog double precision dlog double precision zlog c c The following constants are dlog(2PI)/2, c half, zero, one, seven c data a1, a2, a3, a4, a5/0.918938533204673d0, 0.000595238095238d0, & 0.000793650793651d0, 0.002777777777778d0, & 0.083333333333333d0/ data half, zero, one, seven/0.5d0, 0.0d0, 1.0d0, 7.0d0/ c zlog(f) = dlog(f) c alogam = zero ifault = 1 if (x .lt. zero) return ifault = 0 y = x f = zero if (y .ge. seven) go to 30 f = y 10 y = y + one if (y .ge. seven) go to 20 f = f*y go to 10 20 f = -zlog(f) 30 z = one/(y*y) alogam = f + (y-half)*zlog(y) - y + a1 + (((-a2*z+a3)*z-a4)*z+a5) & /y return end c c*********************************************************************** c EXACT.DAT c*********************************************************************** c 2 8 0.0 0.0 0.0 22 13 5 4 5 3 2 1 7 1 4 3 1 2 3 4 2 10 0.0 0.0 0.0 20 3 6 4 7 6 6 2 2 2 8 8 4 5 2 1 0 2 1 1 3 5 0.0 0.0 0.0 24 7 3 8 1 9 5 5 0 3 2 0 2 0 1 3 7 0.0 0.0 0.0 1 8 5 4 4 2 2 5 3 3 4 3 1 0 10 1 4 0 0 0 0 4 5 0.0 0.0 0.0 9 1 3 1 3 3 3 2 1 3 8 2 4 2 0 4 9 2 0 0 2 5 0.0 0.0 0.0 2 3 4 8 9 0 0 11 10 11 3 3 0.0 0.0 0.0 2 4 6 7 6 1 5 0 0 4 6 0.0 0.0 0.0 2 0 1 2 6 5 1 3 1 1 1 2 1 0 3 1 0 0 1 2 1 2 0 0 5 6 0.0 0.0 0.0 1 2 2 1 1 0 2 0 0 2 3 0 0 1 1 1 2 7 1 1 2 0 0 0 0 1 1 1 1 0 3 9 0.0 0.0 0.0 1 1 1 0 0 0 1 2 4 4 4 4 5 5 5 6 5 0 1 1 1 0 0 0 1 2 4 3 8 0.0 0.0 0.0 20 20 2 2 0 0 0 1 10 10 1 1 3 3 1 0 20 20 2 2 0 0 0 1 3 6 0.0 0.0 0.0 12 6 12 1 1 0 5 12 4 4 0 1 5 12 10 1 1 0 0 0 0.0 0.0 0.0 c c*********************************************************************** c COCHRAN.DAT c*********************************************************************** c 2 8 5.0 80.0 1.0 22 13 5 4 5 3 2 1 7 1 4 3 1 2 3 4 2 10 5.0 80.0 1.0 20 3 6 4 7 6 6 2 2 2 8 8 4 5 2 1 0 2 1 1 3 5 5.0 80.0 1.0 24 7 3 8 1 9 5 5 0 3 2 0 2 0 1 3 7 5.0 80.0 1.0 1 8 5 4 4 2 2 5 3 3 4 3 1 0 10 1 4 0 0 0 0 4 5 5.0 80.0 1.0 9 1 3 1 3 3 3 2 1 3 8 2 4 2 0 4 9 2 0 0 2 5 5.0 80.0 1.0 2 3 4 8 9 0 0 11 10 11 3 3 5.0 80.0 1.0 2 4 6 7 6 1 5 0 0 4 6 5.0 80.0 1.0 2 0 1 2 6 5 1 3 1 1 1 2 1 0 3 1 0 0 1 2 1 2 0 0 5 6 5.0 80.0 1.0 1 2 2 1 1 0 2 0 0 2 3 0 0 1 1 1 2 7 1 1 2 0 0 0 0 1 1 1 1 0 3 9 5.0 80.0 1.0 1 1 1 0 0 0 1 2 4 4 4 4 5 5 5 6 5 0 1 1 1 0 0 0 1 2 4 3 8 5.0 80.0 1.0 20 20 2 2 0 0 0 1 10 10 1 1 3 3 1 0 20 20 2 2 0 0 0 1 3 6 5.0 80.0 1.0 12 6 12 1 1 0 5 12 4 4 0 1 5 12 10 1 1 0 0 0 0.0 0.0 0.0