C MVELPS.FOR Revised 18JUN2001 C MVELPS.FOR 27MAY99 AUTHOR: PAUL N. SOMERVILLE C MVELPS.FOR is a Fortran 90 program to evaluate multivariate normal and multi- C variate-t integrals over a general ellipsoid. Potential known applications C are in reliability for the computation of tolerance factors, and for the C computation of probabilities for linear combinations of (central and non- C central) chisquares and F distributions. C A special case is the computation of cumulative probabilities for noncentral C F and chisquare distributions (for which a number of programs are available). C The program uses a combination of Monte Carlo and quadrature. The procedure C is described in "Numerical Computation of Multivariate Normal and Multi- C variate-t Probabilities over Ellipsoidal Regions" which has been submitted C for publication. A previous publication "Numerical Computation of Multi- C variate Normal and Multivariate-t probabilities over Convex Regions" appears C in the Journal of Computational and Graphical Statistics, vol 7, 1998, pages C 529-524. The present program extends the results to general ellipsoidal reg- C ions in k dimensions. C C Let x be a random vector variable having multivariate normal distribution C with vector mean rvmu and variance covariance matrix vcc = S v where S is C a known k x k positive definite matrix and v is a scalar (commonly called C sigma squared). If v is estimated by s2 with ndenom degrees of freedom, C then x has the multivariate-t distribution. C C The program calculates the probability content of an ellipsoid C (x-emx)' vxinv (x-emx) = ellipg C where emx is the center of the ellipsoid and vxinv is a symmetric k x k C matrix. C C The program transforms (Cholesky transformation) the random vector x C into the product of k uncorrelated normal random variables with unit C variances or a spherically symmetric multivariate-t distribution. C C The program makes "irep" estimates, each of which uses "mocar" random C directions. These are used to obtain an overall estimate and a standard C error of the estimate. A total of 10000 random directions will usually be C sufficient for three decimal accuracy (e.g. irep = 10; mocar = 1000). C C The program can use either of two methods, "binning" or "nonbinning". C The "binning" method is always faster for the multivariate normal case. C The "binning" method is faster for the multivariate-t case when the C number of random directions is large (say ge. 10**6). Estimates and C standard errors for the two methods are nearly identical. C To calculate noncentral F or noncentral chisquare probabilities a shortcut C procedure is available. C C TERMS USED IN THE PROGRAM C k dimension of random vector x C seed seed for random generator C ndenom degrees of freedom for s2, the estimator of sigma squared C (if sigma squared is assumed known, use ndenom = -1). C mocar number of random directions used for each estimate C irep number of independent estimates made (program output is C average of these estimates) C ellipg RHS (scalar) of ellipsoid equation C meth method used to obtain output C for binning method meth = 1 C for nonbinning method meth = 0 C for noncentral "shortcut" meth = -1 C vcc variance covariance matrix of random vector x C rvmu 1 x k mean of random vector x C vxinv k x k matrix of ellipsoid equation (symmetric) C emx 1 x k vector of mean of ellipsoid C znoncent noncentrality parameter (used only for meth = -1) C (Note: "distance" and not "distance squared" is used) C f0 as in Prob[F < f0] for noncentral 'shortcut' method C C RUNNING THE PROGRAM C The program requires that two files exist before running the program. The C file ELIP.IN contains the instructions for the program, and ELIP.OUT is C used for the program output. ELIP.IN can contain several sets of instruct- C ions. After the last set of instructions the program expects a line of 4 or C more -1's. C C CONTENTS OF ELIP.IN C Except for the "shortcut method" used for calculating noncentral chisquare C or noncentral F, (described later), the file ELIP.IN requires six lines of C input. A "line" may actually be a matrix or a vector. The six "lines" C are: C LINE 1 C k seed ndenom mocar irep meth C LINE 2 C ellipg C LINE 3 C vcc C LINE 4 C rvmu C LINE 5 C vxinv C LINE 6 C emx C C LINES 1 and 2: C The elements in lines 1 and 2 have already been defined. If sigma C squared is assumed known, ndenom is given the value of "-1". C LINES 3 and 5: C vcc and vxinv are symmetric k x k matrices. It is permissible to C input only the lower triangular portion of these matrices. Each row C occupies a separate line. (If vcc or vxinv is a diagonal matrix with all C elements equal to "d", a single line containing only "-d" may be used.) C LINES 4 and 6 C the elements of the vectors rmvu and emx make up the lines 4 and 6 C respectively. C C SHORTCUT METHOD FOR NON-CENTRAL PROBABILITIES C The "shortcut method" for calculating noncentral F and noncentral C chisquare probabilities requires two lines only: C LINE 1: k seed ndenom mocar irep meth C LINE 2: znoncent f0 C Lines 1 is the same as before except that now "meth" is set to "-1". C The noncentrality parameter znoncent is defined as a distance (not C distance squared as it sometimes defined!). C f0 is the upper limit for the noncentral F or chisquare as in C Pr[F.le.f0]. C EXAMPLES C (All examples have been compiled with a Lahey FORTRAN 95 compiler and run C using an AMD 800 processor.) C Suppose we have a multivariate-t distribution with k = 3, all correlations C equal to .5, mean vector (0, 0, 0) and degrees of freedom 20. Let the C ellipsoid be given by: C (x - emx)' vxinv (x - emx) = 1 C where C emx = ( 1 .5 .25) C C .25 0 0 C vxinv = 0 .2 0 C 0 0 .1 C C We wish to calculate the required probability first using the binning C method and then using the nonbinning method. C Then ELIP.IN may be given by: C C 3 357 20 10000 10 1 C 1 C 1 C .5 1 C .5 .5 1 C 0 0 0 C .25 C 0 .2 C 0 0 .1 C .1 .5 .25 C C 3 357 20 10000 10 0 C 1 C 1 C .5 1 C .5 .5 1 C 0 0 0 C .25 C 0 .2 C 0 0 .1 C .1 .5 .25 C -1 -1 -1 -1 -1 C Note that the line of negative 1's indicates there are no further sets of C instructions to follow. C The output file ELIP.OUT becomes: C date and time 20010619160723.220 C Binning method used C elapsed time in seconds is 0.390 C no of pops is 3 seed is 357 C df for var est is 20 C value of integral is calculated 10 times C mean value is 0.794763803 C standard error of mean is 4.63318138E-04 C -------------------- C date and time 20010619160723.610 C elapsed time in seconds is 0.930 C no of pops is 3 seed is 357 C df for var est is 20 C value of integral is calculated 10 times C mean value is 0.794791400 C standard error of mean is 4.84645134E-04 C -------------------- C EXAMPLE 2 C Find the probability content of the ellipsoid C (x - emx)'vxinv(x - emx)=1 C where x is MVN(0,I), k=4, emx=(1,1,1,1)' and vxinv=.2I. C We make 10 estimates of the probability, each with 1000 random C directions. We use the binning method. C ELIP.IN may be given by C 4 297 -1 1000 10 1 C 1 C -1 C 0 0 0 0 C -.2 C 1 1 1 1 C -1 -1 -1 -1 -1 -1 C The output file ELIP.OUT becomes: C date and time 20010619160724.540 C Binning method used C elapsed time in seconds is 0.000 C no of pops is 4 seed is 297 C variance assumed known C no of random dir. is 2000.00000 C value of integral is calculated 10 times C mean value is 0.308498234 C standard error of mean is 6.20856299E-04 C -------------------- C EXAMPLE 3 C Find Pr[F<3] where F is the noncentral F distribution with degrees C of freedom 4 and 25 for the numerator and denominator respectively, C and noncentrality parameter 1.2. Make 12 estimates with 1200 random C directions for each. C ELIP.IN may be given by C 4 9133 25 1200 12 -1 C 1.2 3 C -1 -1 -1 -1 -1 -1 C The output file ELIP.OUT becomes: C Noncentrality parameter is 1.20000005 C Upper limit is 3.00000000 C date and time 20010619160724.540 C elapsed time in seconds is 0.220 C no of pops is 4 seed is 9133 C df for var est is 25 C value of integral is calculated 12 times C mean value is 0.905823290 C standard error of mean is 9.21626452E-06 C -------------------- C EXAMPLE 4 C Find Pr[CHI<2.5] where CHI is the noncentral chisquare distribution C with 5 degrees of freedom and noncentrality parameter 1.6. Use a C single estimate with 10000 random directions. C ELIP.IN may be given by C 5 9751 -1 10000 1 -1 C 1.6 2.5 C -1 -1 -1 -1 -1 -1 C The output file ELIP.OUT becomes: C Noncentrality parameter is 1.60000002 C Upper limit is 2.50000000 C date and time 20010619160724.760 C elapsed time in seconds is 0.170 C no of pops is 5 seed is 9751 C variance assumed known C no of random dir. is 100000.000 C value of integral is calculated 1 times C mean value is 9.29978713E-02 C -------------------- C A number of subroutines used in the program have been adapted from C Numerical Methods, by Press et al (Cambridge University Press). C Marsaglia's subroutine rnor is used to generate random normal variates. module comm character(len=28)elapse character(len=8)date character(len=10)time real,dimension(:,:),allocatable::vcc,cmc,vxinv real,dimension(:,:),allocatable::eigvec real,dimension(:,:),allocatable::eigvinv,vw,tinv,vx real,dimension(:),allocatable::dc,wl,pl,prop real, dimension(:),allocatable::qans real,dimension(:),allocatable::emx,emu,eig,rvmu,emeig real,dimension(:),allocatable::preeig,preemeig,preemu integer,dimension(:),allocatable::ibin real sb,zmd2,zmnd2,zmdzn,znd2,fconst2,zm,zn integer k,ndenom,lines,ie,ib,mocar,irep,mm,istart,iend real zero,zone,two,pi,zoneplus,sx,ssx,zrep,sd,zmean real dend integer iran(97),iy,irr,jseed,nz,nw,nran,seed real csum,boundl,boundu,points real cc,e1,t2e,four_e,ellipg,f0 integer ncount1,ncount2,count_rate,count,ier,count_max,iseed parameter (mr=714025,ia=1366,ic=150889,rm=1./mr) parameter (nq=20) real znoncent integer meth real,dimension(:),allocatable::xg,wg,vo parameter(ng=20) end module comm use comm C uses portable random number generator C vc input var-covar matrix C qans vector of q values (irep of them) real se,seconds,poff2 elapse=' elapsed time in seconds is ' Open(8,file='elip.out',access='sequential',form='formatted', $recl=80,status='old') Open(9,file='elip.in',access='sequential',form='formatted', $recl=80,status='old') 21 read(9,*)k,seed,ndenom,mocar,irep,meth if(ndenom.gt.70)then write(8,*)'degrees of freedom exceeds 70. STOP' write(6,*)'degrees of freedom exceeds 70. STOP' go to 997 end if write(6,*)k,seed,ndenom,mocar,irep,meth if((k.lt.0).and.(seed.lt.0).and.(ndenom.lt.0).and.(mocar.lt.0) &)go to 997 if(meth.lt.0)then read(9,*)znoncent,f0 C If k=1 and meth=-1, Monte Carlo not needed. Method is "exact". if(k.eq.1)then mocar=1 irep =1 end if else read(9,*)ellipg end if C set miscellaneous constants zero=0. zone=1. two=2. pi=3.14159265359 iseed=abs(seed) jseed =iseed C iy set to a fixed value for random # generator iy=0 C fill random number bin do j=1,97 iseed=mod(ia*iseed+ic,mr) iran(j)=iseed enddo nz=iran(27) nw=iran(57) nran=iran(97) C irr is seed for racos irr=iseed zm=k zn=ndenom zmd2=zm/two zmnd2=(zm+zn)/two zmdzn=zm/zn znd2=zn/two if(zn.gt.zero)glnu=exp(gammln(zn/two)) 90 sx=0. ssx=0. zrep=irep allocate(qans(irep)) qans=zero if(meth.eq.-1)then C Shortcut for noncentral f and chisq write(8,*)'Noncentrality parameter is',znoncent write(8,*)'Upper limit is ',f0 call date_and_time(date,time) call system_clock(count,count_rate,count_max) ncount1=count write(8,*)'date and time ',date,time do njj=1,irep if(ndenom.gt.0)then call znonf else call znonchi end if qans(njj)=zmean cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm end do deallocate(qans) C End of noncentral calc. Go to rest of printout. go to 910 end if allocate(cmc(k,k),vcc(k,k)) vcc=zero C Input variance covariance matrix (lower triangular portion only) C Input one row at a time C This is varcov for random variable x, read(9,*)vcc(1,1) C Check to see if diagonal matrix if(vcc(1,1).lt.zero)then vcc(1,1)=abs(vcc(1,1)) do i=1,k vcc(i,i)=vcc(1,1) end do else do i=2,k read(9,*) (vcc(i,j),j=1,i) end do do i=1,k-1 do j=i+1,k vcc(i,j)=vcc(j,i) end do end do end if C Input mean vector of random variable allocate(rvmu(k)) read(9,*)(rvmu(i),i=1,k) C***************************************************** call chol C***************************************************** deallocate(vcc) C vcc is now decomposed into (cmc) times (cmc)' C************************************************ C input parameters for the ellipsoidal region of integration C ellipsoid is (x-m)'(inverse of vx)(x-m) allocate(vxinv(k,k)) vxinv=zero C input lower triangular portion of inverse of vx read(9,*)vxinv(1,1) if(vxinv(1,1).le.zero)then vxinv(1,1)=abs(vxinv(1,1)) do i=1,k vxinv(i,i)=vxinv(1,1) end do else do i=2,k read(9,*) (vxinv(i,j),j=1,i) end do do i=1,k-1 do j=i+1,k vxinv(i,j)=vxinv(j,i) end do end do do i=1,k-1 do j=i+1,k vxinv(i,j)=vxinv(j,i) end do end do end if C Input coordinates for ellipsoid center allocate(emx(k)) read(9,*)(emx(i),i=1,k) C make RHS of ellipsoid eqn equal to 1 vxinv=vxinv/ellipg allocate(vx(k,k)) C obtain vx call inverse(vxinv,vx,k) deallocate(vxinv) C Make translation so that random variable x has mean at origin emx=emx-rvmu deallocate (rvmu) allocate(tinv(k,k)) C obtain inverse tinv of cmc call inverse(cmc,tinv,k) deallocate(cmc) C we use transformation w=tinv x to make w into MVN(0,I) C transform vx to vw (w cords) allocate(vw(k,k)) vw=matmul(matmul(tinv,vx),transpose(tinv)) deallocate (vx) allocate(eigvec(k,k)) allocate(eig(k)) C get eigenvalues eig and normalized eigenvectors eigvec call eigprog(vw,eig,eigvec,k,ifault) C write(6,*)'eigenvalues are ',eig C write(6,*)'eigenvectors are ',eigvec C rotate axes so that they parallel ellipse axes C (could have used u=eigvec w) C get ellipse center emu in u axes allocate(eigvinv(k,k)) call inverse(eigvec,eigvinv,k) deallocate(vw) allocate(emu(k)) emu=matmul(matmul(eigvinv,tinv),emx) deallocate(emx,tinv,eigvec,eigvinv) C If we have a spheroid with center at the origin, call sphere 360 if((maxval(eig).eq.minval(eig)).and.(sum(emu).eq.zero))then call sphere(sqrt(eig(1)),poff) write(8,*)'region is a spheroid with center at origin' write(8,*)'value of integral is',poff write(8,*)'--------------------' write(8,*) deallocate(eig,emu,qans) go to 950 end if C Note eqn of ellipsoid is now (u-emu)'(inverse(VV)(u-emu)=1. C where VV has eigenvalues on diag, zeros elsewhere C and u is MVN(0,I) or spherically symmetric multivariate t. C write(6,*)'ellipse center in u coords is ',emu allocate(preeig(k),preemeig(k),emeig(k),preemu(k)) call date_and_time(date,time) call system_clock(count,count_rate,count_max) ncount1=count write(8,*)'date and time ',date,time emeig=emu/eig cac=sum(emu*emeig) preemu=emu preeig=eig preemeig=emeig precac=cac allocate(dc(k)) if(meth.eq.1)then C Binning method begins here C call bound(eig,emu,k,boundl,boundu) C Note eig and emu may come back reordered! C If origin is .lt. .02 from origin, change from "binning" C to "standard' method! C if(boundl.lt..02)then C meth=0 C write(8,*)'Origin is .lt. .02 from boundary' C write(8,*)'Method changed to standard' C go to 480 C end if write(8,*)'Binning method used' if(ndenom.gt.0)then allocate(xg(ng),wg(ng),vo(ng)) do njj=1,irep C sigmasquare is estimated by s2 alf=(zn-two)/two call gaulag(xg,wg,ng,alf) jocar=mocar/ng+1 allocate(pl(nq),wl(nq),prop(nq),ibin(nq)) do ig=1,ng C Get s2 corresponding to xg(ig) s2=xg(ig)*two/zn eig=preeig*s2 emu=preemu call bound(eig,emu,k,boundl,boundu) C Since eig and emu may come back reordered, redo eig! eig=preeig*s2 emu=preemu emeig=preemeig/s2 cac=precac/s2 cac=cac-zone ibin=0 if(cac.le.zero)then C Origin is in ellipsoid xxx=zm/two qsphere=boundl*boundl/two if(qsphere.ge.xxx)then poff=gammp(xxx,qsphere) else poff=zone-gammq(xxx,qsphere) end if call binning4(jocar,cac) vo(ig)=content() vo(ig)=vo(ig)+poff else C Origin is not in ellipsoid call binning5(jocar,cac) vo(ig)=content() end if C******************************************** end do deallocate(pl,wl,prop,ibin) qans(njj)=sum(wg*vo) qans(njj)=qans(njj)/glnu cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm end do deallocate(xg,wg,vo) end if C Insert sigma known , binning here if(ndenom.lt.0)then call bound(eig,emu,k,boundl,boundu) eig=preeig emu=preemu emeig=preemeig cac=cac-zone allocate(pl(nq),wl(nq),prop(nq),ibin(nq)) do njj=1,irep ibin=0 jocar=mocar if(cac.le.zero)then C Origin is in ellipsoid call sphere(boundl,poff) call gauleg(boundl,boundu,pl,wl,nq) do iii=1,jocar call racos aaa=sum(dc*dc/eig) bab=abs(sum(dc*emeig)) b2ac=sqrt(bab*bab-aaa*cac) dis1=abs(bab+b2ac)/aaa call elipbin1(dis1) dis2=abs(bab-b2ac)/aaa call elipbin1(dis2) end do call binprop2(jocar,cac) qans(njj)=content() qans(njj)=qans(njj)+poff else C Origin is not in ellipsoid call binning5(jocar,cac) qans(njj)=content() end if C************************************* cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm end do deallocate(pl,wl,prop,ibin) end if go to 901 end if C Binning method ends here C Non-binning method begins here 480 if(ndenom.gt.0)allocate(xg(ng),wg(ng),vo(ng)) do 900 njj=1,irep if(ndenom.gt.0) then C sigmasquare is estimated by s2 alf=(zn-two)/two call gaulag(xg,wg,ng,alf) jocar=mocar/ng+1 do ig=1,ng C Get s2 corresponding to xg(ig) s2=xg(ig)*two/zn eig=preeig*s2 emeig=preemeig/s2 cac=precac/s2 cac=cac-zone summ=zero do iii=1,jocar call racos aaa=sum(dc*dc/eig) bab=sum(dc*emeig) b2ac=bab*bab-aaa*cac if(b2ac.ge.zero)then C Line goes thru ellipsoid b2ac=sqrt(b2ac) call gamtrib(poff1,poff2,aaa,bab,b2ac) if(cac.le.zero)then C Center is in ellipsoid summ=summ +poff1+poff2 else C Center is not in ellipsoid summ=summ+abs(poff1-poff2) end if end if end do vo(ig)=summ/jocar/two end do points=jocar*ng qans(njj)=sum(wg*vo) qans(njj)=qans(njj)/glnu end if C Still non-binning method here if(ndenom.lt.0) then C Variance is unknown summ=zero cac=precac-zone do i=1,mocar call racos aaa=sum(dc*dc/eig) bab=sum(dc*emeig) b2ac=bab*bab-aaa*cac if(b2ac.ge.zero)then C Line goes thru ellipsoid b2ac=sqrt(b2ac) call gamtrib(poff1,poff2,aaa,bab,b2ac) if(cac.le.zero)then C Center is in ellipsoid summ=summ+poff1+poff2 else C Center is not in ellipsoid summ=summ+abs(poff1-poff2) end if end if 610 end do points=mocar*2 qans(njj)=summ/points end if cc=qans(1) cm=qans(njj)-cc sx=sx+cm ssx=ssx+cm*cm 900 continue if(ndenom.gt.0)deallocate(xg,wg,vo) 901 deallocate(eig,emeig,emu,dc) deallocate(preeig,preemeig,preemu) deallocate(qans) C End of non-binning method 910 call date_and_time(date,time) call system_clock(count) ncount2=count seconds=ncount2-ncount1 if(seconds.lt.0)seconds=seconds+count_max seconds=seconds/count_rate write(8,8989)elapse,seconds 8989 format(a28,f10.3) sd=9999. se=9999. if(zrep.eq.1) go to 920 sd=(ssx-sx*sx/zrep)/(zrep-1.) if(sd.le.zero)sd=1.E-20 se=sqrt(sd/zrep) sd=sqrt(sd) sd=sd/two se=se/two 920 zmean=sx/zrep+cc write(8,*)'no of pops is ',k,' seed is',jseed if(ndenom.gt.0)then write(8,*)' df for var est is',ndenom else write(8,*)' variance assumed known' write(8,*)'no of random dir. is',points end if write(8,*)'value of integral is calculated', irep,' times' write(8,*)'mean value is ',zmean if(zrep.eq.zone)go to 940 write(8,*)'standard error of mean is ',se 940 write(8,*)'--------------------' 9753 format(1x,f3.0,f4.0,i3,f9.6,f9.7,f10.3) 950 continue C WE ARE NOW FINISHED AND NEED TO GO BACK for new instructions GO TO 21 997 CLOSE(8) CLOSE(9) 998 stop 8000 format(10(f6.4,1x)) end C**************************************** subroutine binning4(jocar,cac) use comm call gauleg(boundl,boundu,pl,wl,nq) do iii=1,jocar call racos aaa=sum(dc*dc/eig) bab=abs(sum(dc*emeig)) b2ac=sqrt(bab*bab-aaa*cac) dis1=abs(bab+b2ac)/aaa call elipbin1(dis1) dis2=abs(bab-b2ac)/aaa call elipbin1(dis2) end do call binprop2(jocar,cac) return end C************************************** subroutine binning5(jocar,cac) use comm call gauleg(boundl,boundu,pl,wl,nq) do iii=1,jocar call racos aaa=sum(dc*dc/eig) bab=abs(sum(dc*emeig)) b2ac=bab*bab-aaa*cac if(b2ac.ge.zero)then C Line goes thru ellipsoid b2ac=sqrt(b2ac) dis1=(bab+b2ac)/aaa call elipbin1(dis1) dis2=(bab-b2ac)/aaa call elipbin2(dis2) end if end do call binprop2(jocar,cac) return end C*************************************** subroutine gamtrib(poff1,poff2,aaa,bab,b2ac) C (Routine obtains elements for probability calc) use comm xxx=zm/two t=(bab+b2ac)/aaa qsphere=t*t/two if(qsphere.ge.xxx)then poff1=gammp(xxx,qsphere) else poff1=zone-gammq(xxx,qsphere) end if t=(bab-b2ac)/aaa qsphere=t*t/two if(qsphere.ge.xxx)then poff2=gammp(xxx,qsphere) else poff2=zone-gammq(xxx,qsphere) end if return end subroutine znonf use comm allocate(dc(k)) allocate(xg(ng),wg(ng),vo(ng)) alf=(zn-two)/two call gaulag(xg,wg,ng,alf) jocar=mocar/ng+1 points=jocar*zm*two do ig=1,ng C Get s2 corresponding to xg(ig) summ=zero s2=xg(ig)*two/zn ta=znoncent**2-zm*s2*f0 do iii=1,jocar call racos do i=1,k tdc=znoncent*dc(i) te=tdc**2-ta if(te.ge.zero)then C Line goes thru ellipsoid te=sqrt(te) call tribcen(poff1,poff2,tdc,te) if(ta.le.zero)then C Ellipsoid contains origin summ=summ+poff1+poff2 else C Ellipsoid does not contain origin summ=summ+abs(poff1-poff2) end if end if end do end do vo(ig)=summ/points end do glnu=exp(gammln(zn/two)) zmean=sum(vo*wg)/glnu deallocate(dc) deallocate(xg,wg,vo) return end subroutine znonchi use comm allocate(dc(k)) summ=zero ta=znoncent**2-f0 points=mocar*zm*two do iii=1,mocar call racos do i=1,k tdc=znoncent*dc(i) te=tdc**2-ta if(te.ge.zero)then C Line goes thru ellipsoid te=sqrt(te) call tribcen(poff1,poff2,tdc,te) if(ta.le.zero)then C Ellipsoid contains origin summ=summ+poff1+poff2 else C Ellipsoid does not contain origin summ=summ+abs(poff1-poff2) end if end if end do end do zmean=summ/points deallocate(dc) return end subroutine tribcen(poff1,poff2,t1,t2) use comm xxx=zm/two t=t1+t2 qsphere=t*t/two if(qsphere.ge.xxx)then poff1=gammp(xxx,qsphere) else poff1=zone-gammq(xxx,qsphere) end if t=t1-t2 qsphere=t*t/two if(qsphere.ge.xxx)then poff2=gammp(xxx,qsphere) else poff2=zone-gammq(xxx,qsphere) end if return end subroutine elipbin1(dist) C Increments distances to proper bin use comm do 40 i=1,nq if(dist.ge.pl(i))go to 40 ibin(i)=ibin(i)+1 go to 50 40 end do 50 return end subroutine elipbin2(dist) C Decrements distances to proper bin use comm do 40 i=1,nq if(dist.ge.pl(i))go to 40 ibin(i)=ibin(i)-1 go to 50 40 end do 50 return end C**************************************** C This routine calculates proportions in the bins subroutine binprop2(kocar,cac) use comm integer i points=two*kocar do i=2,nq ibin(i)=ibin(i-1)+ibin(i) end do prop=ibin prop=-prop/points if(cac.le.zero)prop=zone+prop return end C*************************************** function content() C Obtains prob content for an ellipsoid corr. to an s2 use comm ce=two**(zmd2-zone)*exp(gln(zmd2)) content=zero do i=1,nq cd=pl(i)**(zm-zone)*exp(-pl(i)**2/two) content=content+cd*wl(i)*prop(i) end do content=content/ce return end C**************************************** Subroutine pf(n1,n2,z,poff,ier) C (This is a routine to obtain probabilities for the F-distribution) ier=0 if(z .gt. 0.0d0) go to 5 poff=0.0d0 ier=1 return 5 if(n1.gt.0 .and. n2.gt.0) go to 10 ier=2 poff=0.0 return 10 continue an1=n1 an2=n2 a=an1*z/(an1*z+an2) a1=1.0d0-a if(a1 .lt. 0.1d-36) a1=0.1d-36 d1=an1*0.5d0 d2=an2*0.5d0 d3=d1+d2-1.0d0 r=0.d0 s1=0.d0 s2=0.d0 del=1.0d0 xm=1.0d0 xk=1.0d0 cc=0.25d0 pi=3.141592653589793d0 n=n2 C NOTE BEGINNING OF MAJOR LOOP 15 continue C TO SEE IF DEGREES OF FREEDOM ARE ODD OR EVEN m=d2 m=2*m if(m .ne. n) go to 30 n=d2-1 C IF DEGREES OF FREEDOM ARE EVEN C N=D.F./2-1 if(n .eq. 0) go to 25 do 20 i=1,n s1=del+s1*r d2=d2-1.0d0 d3=d3-1.0d0 tem=a1/d2 r=d3*tem s2=(r+tem)*s2 20 continue 25 s1=del+s1*r del=0.d0 t=-1.0d0 d3=-1.0d0 s2=a*s2 cc=cc+0.5d0 go to 45 C IF DEGREES OF FREEDOM ARE ODD C N=(D.F.-1)/2 30 n=d2 C IF DEGREES OF FREEDOM EQUAL 1. C DO NO EXIT LOOP. if(n .eq. 0)go to 40 do 35 i=1,n s1=del+s1*r d2=d2-1.0d0 d3=d3-1.0d0 tem=a1/d2 r=d3*tem s2=(r+tem)*s2 35 continue 40 s1=xk*s1 s2=xk*s2 art=sqrt(a1) xm=xm*art t=(xm-art)/a1 d3=-0.5d0 xk=2.0d0/pi cc=cc*2.0d0 45 if(cc .gt. 0.875d0) go to 50 d2=d1 d3=d2+d3 s2=s1 s1=0.0d0 a1=a if(a1 .lt. 0.1d-36) a1=0.1d-36 n=n1 go to 15 50 if(cc .lt. 1.125d0) del=4.0d0/pi*atan(t) poff=xm*(s2-s1)-del if(0.0d0 .le. poff .and. 1.0d0 .ge. poff) return if(poff .lt. 0.0d0) poff=0.0d0 if(poff .gt. 1.0d0) poff=1.0d0 ier=3 return end C**************************************** C Subroutine to generate direction cosines of a random direction in k C dimensions subroutine racos use comm real sum 10 sum=0. do 20 i=1,k dc(i)=rnor(nz,nw) sum=sum+dc(i)*dc(i) 20 continue sum=sqrt(sum) do i=1,k dc(i)=dc(i)/sum end do return end C**************************************** function gln(a) C This function gives the logarithm of the gamma function when the C argument is a positive integer or an integer divided by 2 implicit real (a-h,o-z) i1=aint(a) c2=0. c1=1. gln=0. if((a-i1).eq.0.) then i1=i1-2 else c1=-.5 gln=alog(sqrt(acos(-1.))) endif do 20 i=1,i1 z=c1+i x=alog(z) gln=gln+x 20 continue return end C******************************************** function beta(a,b) C This function gives the ratio of gamma(a+b) to gamma(a)*gamma(b) beta=gln(a+b)-gln(a)-gln(b) beta= exp(beta) return end C************************************************* subroutine chol C decomposes pos def matrix vcc into cmc times cmc transpose use comm real,dimension(k)::dd integer i,ii,j cmc=0. do 20 i=1,k cmc(i,i)=1. 20 continue do 100 j=1,k dd(j)=vcc(j,j) do 50 ii=1,j-1 dd(j)=dd(j)-cmc(j,ii)*cmc(j,ii)*dd(ii) 50 continue do 80 i=j+1,k cmc(i,j)=vcc(i,j) do 60 ii=1,j-1 cmc(i,j)=cmc(i,j)-cmc(j,ii)*cmc(i,ii)*dd(ii) 60 continue cmc(i,j)=cmc(i,j)/dd(j) 80 continue 100 continue do 120 j=1,k a=sqrt(dd(j)) do 110 i=j,k cmc(i,j)=a*cmc(i,j) 110 continue 120 continue C write(8,*)'cholesky decomposition' C write(8,*)'written from chol' C do i=1,k C write(8,*)(cmc(i,j),j=1,k) C end do return end C***************************************************************** C gammp, gammq,gser,gcf and gammln are from "Numerical Methods". FUNCTION GAMMP(A,X) IF(X.LT.0..OR.A.LE.0.)write(8,*) 'input error in GAMMP' IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMP=GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMP=1.-GAMMCF ENDIF RETURN END FUNCTION GAMMQ(A,X) IF(X.LT.0..OR.A.LE.0.)write(8,*)'input error in GAMMQ' IF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMQ=1.-GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMQ=GAMMCF ENDIF RETURN END SUBROUTINE GSER(GAMSER,A,X,GLN) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) IF(X.LE.0.)THEN IF(X.LT.0.)write(8,*)'input error in GSER' GAMSER=0. RETURN ENDIF AP=A SUM=1./A DEL=SUM DO 11 N=1,ITMAX AP=AP+1. DEL=DEL*X/AP SUM=SUM+DEL IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1 11 CONTINUE write(8,*)'A too large, ITMAX too small in GSER' 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) RETURN END SUBROUTINE GCF(GAMMCF,A,X,GLN) PARAMETER (ITMAX=100,EPS=3.E-7) GLN=GAMMLN(A) GOLD=0. A0=1. A1=X B0=0. B1=1. FAC=1. DO 11 N=1,ITMAX AN=FLOAT(N) ANA=AN-A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF=AN*FAC A1=X*A0+ANF*A1 B1=X*B0+ANF*B1 IF(A1.NE.0.)THEN FAC=1./A1 G=B1*FAC IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1 GOLD=G ENDIF 11 CONTINUE write(8,*)'A too large, ITMAX too small in GCF' 1 GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G RETURN END FUNCTION GAMMLN(XX) REAL COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173,-86.50532033,24.01409822, * -1.231739516,.120858003E-2,-.536382E-5,2.50662827465/ DATA HALF,ONE,FPF/0.5,1.0,5.5/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+LOG(STP*SER) RETURN END Subroutine gauleg(x1,x2,x,w,n) C From "Numerical Methods" p 125. Gives constants for Gauss-Legendre. C Original range is x1 to x2. Abscissas and weights are arrays x and w C of length N. implicit double precision(a-h,o-z) real x1,x2,x(n),w(n) parameter (eps=3.d-14) m=(n+1)/2 xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) do 12 i=1,m z=cos(3.141592654d0*(i-.25d0)/(n+.5d0)) 1 continue p1=1.0d0 p2=0.0d0 do 11 j=1,n p3=p2 p2=p1 p1=((2.0d0*j-1.0d0)*z*p2-(j-1.d0)*p3)/j 11 continue pp=n*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp if(abs(z-z1).gt.eps) go to 1 x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.0d0*xl/((1.d0-z*z)*pp*pp) w(n+1-i)=w(i) 12 continue return end module inv integer,dimension(:),allocatable::indx integer joj end module inv subroutine inverse(vcc,vinv,k) C vinv is calculated as the inverse of vcc C the matrix need not be symmetric use inv real vcc(k,k),vinv(k,k) allocate(indx(k)) vinv=0. do i=1,k vinv(i,i)=1. end do call ludcmp(vcc,k) do joj=1,k call lubksb(vcc,vinv,k) end do deallocate(indx) return end subroutine lubksb(vcc,vinv,k) use inv real vcc(k,k),vinv(k,k) INTEGER ii,ll REAL sum ii=0 do 12 i=1,k ll=indx(i) sum=vinv(ll,joj) vinv(ll,joj)=vinv(i,joj) if (ii.ne.0)then do 11 jj=ii,i-1 sum=sum-vcc(i,jj)*vinv(jj,joj) 11 continue else if (sum.ne.0.) then ii=i endif vinv(i,joj)=sum 12 continue do 14 i=k,1,-1 sum=vinv(i,joj) do 13 jj=i+1,k sum=sum-vcc(i,jj)*vinv(jj,joj) 13 continue vinv(i,joj)=sum/vcc(i,i) 14 continue return END SUBROUTINE ludcmp(vcc,k) use inv real vcc(k,k) real,dimension(:),allocatable::vv INTEGER NMAX REAL d,TINY PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER imax REAL aamax,dum,sum allocate(vv(nmax)) C d is used only for calculating sign of determinant d=1. do 12 i=1,k aamax=0. do 11 joj=1,k if (abs(vcc(i,joj)).gt.aamax) aamax=abs(vcc(i,joj)) 11 continue if (aamax.eq.0.)write(6,*)'singular matrix in ludcmp' vv(i)=1./aamax 12 continue do 19 joj=1,k do 14 i=1,joj-1 sum=vcc(i,joj) do 13 kk=1,i-1 sum=sum-vcc(i,kk)*vcc(kk,joj) 13 continue vcc(i,joj)=sum 14 continue aamax=0. do 16 i=joj,k sum=vcc(i,joj) do 15 kk=1,joj-1 sum=sum-vcc(i,kk)*vcc(kk,joj) 15 continue vcc(i,joj)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (joj.ne.imax)then do 17 kk=1,k dum=vcc(imax,kk) vcc(imax,kk)=vcc(joj,kk) vcc(joj,kk)=dum 17 continue d=-d vv(imax)=vv(joj) endif indx(joj)=imax if(vcc(joj,joj).eq.0.)vcc(joj,joj)=TINY if(joj.ne.k)then dum=1./vcc(joj,joj) do 18 i=joj+1,k vcc(i,joj)=vcc(i,joj)*dum 18 continue endif 19 continue deallocate(vv) return end ************************************************ module meig real,dimension(:,:),allocatable::zzz real,dimension(:),allocatable::eee real,parameter::eps=1.e-14 real,parameter::eta=1.e-37 real,parameter::precis=1.e-14 real zone,zero,two,tol integer,parameter::mits=30 end module meig subroutine eigprog(vcc,eig,eigvec,k,ifault) use meig real vcc(k,k),eig(k),eigvec(k,k) allocate(zzz(k,k),eee(k)) zone=1. zero=0. two=2. call tdiag(vcc,eig,k,ifault) call lrvt(eig,k,ifault) eigvec=zzz deallocate(zzz,eee) return end SUBROUTINE TDIAG(vcc,eig,k,ifault) C C Algorithm as 60.1 appl.statist. (1973) vol.22 no.2 C C reduces real symmetric matrix to tridiagonal form C C tol is a machine dependent constant , tol = eta/eps , where C eta is the smallest positive number representable in the C computer and eps is the smallest positive number for which C 1+eps.ne.1. C C k is input order of symmetric matrix vcc C vcc Real array(k,k). Input. Only the lower triangle is used C eig Real array (k).Output.The diagonal elements of the tridiagonal matrix. C eee Real array (k).Output.Element2 2 to k give the off diagonal elements C of the tridiagonal matrix. Element 1 is set to zero. C zzz Real array (k,k).Output. Product of the Householder transformation. use meig real vcc(k,k),eig(k) REAL F,G,H,HH INTEGER L,J,I,I1,J1,KK TOL=ETA/EPS IFAULT=1 IF (K.LE.1) RETURN IFAULT=0 DO I=1,K DO J=1,I ZZZ(I,J)=vcc(I,J) END DO END DO I=K DO I1=2,K L=I-2 F=ZZZ(I,I-1) G=ZERO IF (L.LT.1) GO TO 30 DO KK=1,L G=G+ZZZ(I,KK)**2 END DO 30 H=G+F*F C C if g is too small for orthogonality to be guaranteed, the C transformation is skipped C IF (G.GT.TOL) GO TO 40 EEE(I)=F EIG(I)=ZERO GO TO 100 40 L=L+1 G=SQRT(H) IF (F.GE.ZERO) G=-G EEE(I)=G H=H-F*G ZZZ(I,I-1)=F-G F=ZERO DO J=1,L ZZZ(J,I)=ZZZ(I,J)/H G=ZERO C C form element of aaa * u C DO KK=1,J G=G+ZZZ(J,KK)*ZZZ(I,KK) END DO IF (J.GE.L) GO TO 70 J1=J+1 DO KK=J1,L G=G+ZZZ(KK,J)*ZZZ(I,KK) END DO C C form element of p C 70 EEE(J)=G/H F=F+G*ZZZ(J,I) END DO C C form k C HH=F/(H+H) C C form reduced vcc C DO J=1,L F=ZZZ(I,J) G=EEE(J)-HH*F EEE(J)=G DO KK=1,J ZZZ(J,KK)=ZZZ(J,KK)-F*EEE(KK)-G*ZZZ(I,KK) END DO END DO EIG(I)=H 100 I=I-1 END DO EIG(1)=ZERO EEE(1)=ZERO C C accumulation of transformation matrices C DO I=1,K L=I-1 IF (EIG(I).EQ.ZERO.OR.L.EQ.0) GO TO 140 DO J=1,L G=ZERO DO KK=1,L G=G+ZZZ(I,KK)*ZZZ(KK,J) END DO DO KK=1,L ZZZ(KK,J)=ZZZ(KK,J)-G*ZZZ(KK,I) END DO END DO 140 EIG(I)=ZZZ(I,I) ZZZ(I,I)=ZONE IF (L.EQ.0) GO TO 160 DO J=1,L ZZZ(I,J)=ZERO ZZZ(J,I)=ZERO END DO 160 CONTINUE END DO RETURN END C C SUBROUTINE LRVT(eig,k,ifault) C C algorithm as 60.2 appl.statist. (1973) vol.22, no.2 C C finds latent roots and vectors of tridiagonal matrix C C k integer input. order of the tridiagonal matrix. C EIG Real (k) Input diagonal elements of the tridiagonal matrix C Output the latent roots C EEE Real (k) Input:Householder transformation output fron TDIAG C if latent vectors of the matrix required. C ZZZ Real (k,k) Output: normalized latent vectors, col by col. C IFAULT Output: fault indicator (1 if more than MITS iterations C 0 otherwise). MITS set to zero,(can be reset!) use meig real eig(k) REAL F,G,H,B,P,R,PR,C,S INTEGER L,J,I,I1,JJ,M,M1,K1,KK IFAULT=2 IF (K.LE.1) RETURN IFAULT=1 K1=K-1 DO I=2,K EEE(I-1)=EEE(I) END DO EEE(K)=ZERO B=ZERO F=ZERO DO 100 L=1,K JJ=0 H=PRECIS*(ABS(EIG(L))+ABS(EEE(L))) IF (B.LT.H) B=H C C look for small sub-diagonal element C DO 20 M1=L,K M=M1 IF (ABS(EEE(M)).LE.B) GO TO 30 20 CONTINUE 30 IF (M.EQ.L) GO TO 90 40 IF (JJ.EQ.MITS) RETURN JJ=JJ+1 C C form shift C P=(EIG(L+1)-EIG(L))/(TWO*EEE(L)) R=SQRT(P*P+ZONE) PR=P+R IF (P.LT.ZERO) PR=P-R H=EIG(L)-EEE(L)/PR DO I=L,K EIG(I)=EIG(I)-H END DO F=F+H C C ql transformation C P=EIG(M) C=ZONE S=ZERO M1=M-1 I=M DO I1=L,M1 J=I I=I-1 G=C*EEE(I) H=C*P IF (ABS(P).GE.ABS(EEE(I))) GO TO 60 C=P/EEE(I) R=SQRT(C*C+ZONE) EEE(J)=S*EEE(I)*R S=ZONE/R C=C/R GO TO 70 60 C=EEE(I)/P R=SQRT(C*C+ZONE) EEE(J)=S*P*R S=C/R C=ZONE/R 70 P=C*EIG(I)-S*G EIG(J)=H+S*(C*G+S*EIG(I)) C C form vector C DO KK=1,K H=ZZZ(KK,J) ZZZ(KK,J)=S*ZZZ(KK,I)+C*H ZZZ(KK,I)=C*ZZZ(KK,I)-S*H END DO END DO EEE(L)=S*P EIG(L)=C*P IF (ABS(EEE(L)).GT.B) GO TO 40 90 EIG(L)=EIG(L)+F 100 CONTINUE C C order latent roots and vectors C C DO 130 I=1,K1 C KK=I C P=EIG(I) C I1=I+1 C DO 110 J=I1,K C IF (EIG(J).LE.P) GO TO 110 C KK=J C P=EIG(J) C110 CONTINUE C IF (KK.EQ.I) GO TO 130 C EIG(KK)=EIG(I) C EIG(I)=P C DO 120 J=1,K C P=ZZZ(J,I) C ZZZ(J,I)=ZZZ(J,KK) C ZZZ(J,KK)=P C120 CONTINUE C130 CONTINUE C IFAULT=0 RETURN END C*********************************** module bou real, dimension(:),allocatable::eig,elm,amda,dist logical,dimension(:),allocatable::ok logical err integer k real elms,elmb,ends,endb real,parameter::tiny=1.e-20 real,parameter::factor=.0000001 real,parameter::tol=.0001 real,parameter::xacc=.00002 end module subroutine bound(teig,telm,kk,boundl,boundu) C Calculates the minimum and maximum distances from the origin C to the boundaries of an ellipsoid. use bou real,dimension(:),allocatable::u real teig(kk),telm(kk) if(dot_product(telm,telm).eq.0.)then boundl=sqrt(minval(teig)) boundu=sqrt(maxval(teig)) go to 999 end if C Above bounds hold when ellipse center is at origin. telm=-telm call sort2(kk,telm,teig) telm=-telm C Find i, the number of elements of telm which are zero do i=kk,1,-1 if(telm(i).ne.0.)then k=i go to 30 end if 20 end do 30 if(k.eq.kk)go to 40 C Find extrema corr. to cases where element of telm .eq.0 allocate(u(k)) tempu=0. templ=999999 C Calculate vector u. C We "go to 33" when we get impossible values for u(i) do j=k+1,kk temp=0. do i=1,k if(teig(i).eq.teig(j))go to 33 u(i)=-teig(j)*telm(i)/(teig(i)-teig(j)) if((u(i).gt.(telm(i)+sqrt(teig(i)))).or. &(u(i).lt.(telm(i)-sqrt(teig(i)))))go to 33 C if((u(i).gt.sqrt(teig(i))+telm(i)).or. C &(u(i).lt.-sqrt(teig(i))+telm(i)))go to 33 temp=temp+(u(i)-telm(i))**2/teig(i) end do if(temp.gt.1.)go to 33 C Value for u(j)**2 temp=teig(j)*(1.-temp) C Value for new extreme distance temp=sqrt(temp+dot_product(u,u)) C write(6,*)'an2 extreme',temp if(temp.gt.tempu)tempu=temp if(temp.lt.templ)templ=temp 33 end do deallocate(u) C Find extrema corr. to cases where telm elements .ne.0 40 allocate(eig(k),elm(k)) eig=teig(1:k) elm=telm(1:k) allocate(amda(2*k),ok(k-1),dist(2*k)) call sort2(k,eig,elm) C write(6,*)'eig',eig C write(6,*)'elm',elm ok=.true. C Find all lamda solutions between -eig(k) and -eig(1) do i=1,k-1 err=.true. ends=-eig(i+1) endb=-eig(i) elms=elm(i+1) elmb=elm(i) if((endb-ends).ge.tiny)then call mybrak(ax,bx,cx) C write(6,*)'ax=',ax,'bx=',bx,'cx=',cx C write(6,*)'OK= ',ok(i) if(err.eqv..false.)ok(i)=.false. else ok(i)=.false. end if fmin=golden(ax,bx,cx,tol,xmin) C write(6,*)'minimum is ',fmin,' at lamda= ',xmin if((ok(i).neqv..false.).and.(fmin.le.0.))then amda(2*i)=rtsafe(xmin,cx,xacc) amda(2*i+1)=rtsafe(ax,xmin,xacc) C write(6,*)'for i=',i,' Roots are ',amda(2*i),amda(2*i+1) else ok(i)=.false. amda(2*i)=0. amda(2*i+1)=0. C write(6,*)'for i= ',i,' no roots' end if end do C Get smallest and largest lamda solutions. ax=-eig(1)+factor*abs(eig(1)) amda(1)=rtsafe(ax,999999.,xacc) cx=-eig(k)-factor*abs(eig(k)) amda(2*k)=rtsafe(-999999.,cx,xacc) C write(6,*)'smallest and largest roots are ',amda(2*k),amda(1) do i=1,k-1 if(ok(i).eqv..true.)then dist(2*i)=(amda(2*i))**2*sum((elm/(eig+amda(2*i)))**2) dist(2*i)=sqrt(dist(2*i)) dist(2*i+1)=(amda(2*i+1))**2*sum((elm/(eig+amda(2*i+1)))**2) dist(2*i+1)=sqrt(dist(2*i+1)) else dist(2*i)=-1 dist(2*i+1)=-1 end if end do dist(1)=(amda(1))**2*sum((elm/(eig+amda(1)))**2) dist(1)=sqrt(dist(1)) dist(2*k)=(amda(2*k))**2*sum((elm/(eig+amda(2*k)))**2) dist(2*k)=sqrt(dist(2*k)) C write(6,*)'distances (max or min) from origin to "ellipsoid" C &in order of increasing lamda values are: - negative distances C &indicate neither max or min at this lamda -' C write(6,*)'lambda solutions are' C write(6,*)amda C write(6,*)'corresponding distances are' C write(6,*)dist C write(6,*)'Zero values for lamda indicate no real solutions' C write(6,*)'Values of -1 for distance indicate no max or min' boundl=9999999. do i=1,k if((dist(i).ne.-1.).and.(dist(i).lt.boundl))boundl=dist(i) end do boundu=maxval(dist) if(k.eq.kk)go to 900 C write(6,*)'templ,tempu',templ,tempu boundl=min(boundl,templ) boundu=max(boundu,tempu) 900 deallocate(amda,ok,dist) deallocate(eig,elm) C 999 write(6,*)'Bounds for distance to the origin are ' C write(6,*)boundl,boundu 999 return end subroutine mybrak(ax,bx,cx) use bou del=factor*ends ax=ends+abs(del) bx=(ends+endb)/2. del=factor*endb cx=endb-abs(del) if((flam(ax).lt.flam(bx)).or.(flam(cx).lt.flam(bx))) &err=.false. return end function flam(alam) use bou flam=sum(eig*(elm/(eig+alam))**2)-1. return end function flamd(alam) use bou flamd=-2.*sum(eig*(elm/(eig+alam))**3) return end SUBROUTINE sort2(n,arr,brr) INTEGER n,M,NSTACK REAL arr(n),brr(n) PARAMETER (M=7,NSTACK=50) INTEGER i,ir,j,jstack,k,l,istack(NSTACK) REAL a,b,temp jstack=0 l=1 ir=n 1 if(ir-l.lt.M)then do 12 j=l+1,ir a=arr(j) b=brr(j) do 11 i=j-1,l,-1 if(arr(i).le.a)goto 2 arr(i+1)=arr(i) brr(i+1)=brr(i) 11 continue i=l-1 2 arr(i+1)=a brr(i+1)=b 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp temp=brr(k) brr(k)=brr(l+1) brr(l+1)=temp if(arr(l).gt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp temp=brr(l) brr(l)=brr(ir) brr(ir)=temp endif if(arr(l+1).gt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp temp=brr(l+1) brr(l+1)=brr(ir) brr(ir)=temp endif if(arr(l).gt.arr(l+1))then temp=arr(l) arr(l)=arr(l+1) arr(l+1)=temp temp=brr(l) brr(l)=brr(l+1) brr(l+1)=temp endif i=l+1 j=ir a=arr(l+1) b=brr(l+1) 3 continue i=i+1 if(arr(i).lt.a)goto 3 4 continue j=j-1 if(arr(j).gt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp temp=brr(i) brr(i)=brr(j) brr(j)=temp goto 3 5 arr(l+1)=arr(j) arr(j)=a brr(l+1)=brr(j) brr(j)=b jstack=jstack+2 if(jstack.gt.NSTACK)then write(6,*)'NSTSCK too small in sort2' end if if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 END FUNCTION golden(ax,bx,cx,tol,xmin) REAL golden,ax,bx,cx,tol,xmin,R,C C tol is the precision of the minimum golden C xkin is the corresponding abscissa PARAMETER (R=.61803399,C=1.-R) REAL f1,f2,x0,x1,x2,x3 x0=ax x3=cx if(abs(cx-bx).gt.abs(bx-ax))then x1=bx x2=bx+C*(cx-bx) else x2=bx x1=bx-C*(bx-ax) endif f1=flam(x1) f2=flam(x2) 1 if(abs(x3-x0).gt.tol*(abs(x1)+abs(x2)))then if(f2.lt.f1)then x0=x1 x1=x2 x2=R*x1+C*x3 f1=f2 f2=flam(x2) else x3=x2 x2=x1 x1=R*x2+C*x0 f2=f1 f1=flam(x1) endif goto 1 endif if(f1.lt.f2)then golden=f1 xmin=x1 else golden=f2 xmin=x2 endif return END FUNCTION rtsafe(x1,x2,xacc) C finds root of function bracketed between x1 and x2 within xacc INTEGER MAXIT REAL rtsafe,x1,x2,xacc PARAMETER (MAXIT=100) INTEGER j REAL df,dx,dxold,f,fh,fl,temp,xh,xl fl=flam(x1) fh=flam(x2) df=flamd(x2) if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))then write(6,*)'root must be bracketed in rtsafe' write(6,*)'rtsafe x1,x2',x1,x2 write(6,*)'fl, fh',fl,fh end if if(fl.eq.0.)then rtsafe=x1 return else if(fh.eq.0.)then rtsafe=x2 return else if(fl.lt.0.)then xl=x1 xh=x2 else xh=x1 xl=x2 endif rtsafe=.5*(x1+x2) dxold=abs(x2-x1) dx=dxold f=flam(rtsafe) df=flamd(rtsafe) do 11 j=1,MAXIT if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0..or. abs(2.* *f).gt.abs(dxold*df) ) then dxold=dx dx=0.5*(xh-xl) rtsafe=xl+dx if(xl.eq.rtsafe)return else dxold=dx dx=f/df temp=rtsafe rtsafe=rtsafe-dx if(temp.eq.rtsafe)return endif if(abs(dx).lt.xacc) return f=flam(rtsafe) df=flamd(rtsafe) if(f.lt.0.) then xl=rtsafe else xh=rtsafe endif 11 continue write(6,*)'WARNING:max iterations exceeded in rtsafe' write(6,*)'Diff of ',dx,' Target was +or- ',xacc return END ********************************** subroutine sphere(q,poff) use comm C Find probability within sphere if(ndenom.gt.0)then qsphere=q*q/zm call pf(k,ndenom,qsphere,poff,ier) else xxx=zm/two qsphere=q*q/two if(qsphere.ge.xxx)then poff=gammp(xxx,qsphere) else poff=zone-gammq(xxx,qsphere) end if end if return end function gamma(a,z,w,jran) integer z,w,jran real a d=a-.3333333 C c=.1111111/sqrt(d) c=.3333333/sqrt(d) 1 x=rnor(z,w) v=1.+c*x if(v .le. 0.) go to 1 v=v*v*v gamma=d*v jran=69069*jran U=.5+jran*.2328306e-9 C if( x.le.2. .and. U.le.1.-.12936*x*x) return if(U.le.1.-.0331*x**4)return if(log(U)-.5*x*x .le. d*(1.-v+log(v))) return go to 1 return end c*** Monty Python Normal function rnor(z,w) integer z,w data a,b,s/1.17741,2.506628,.8857913/ integer,parameter::m1=18000,m2=65535 z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) c*** rnor uniform, -b to b, via b*.5^31*IVNI rnor=1.16724e-9*(ishft(z,16)+iand(w,m2)) if(abs(rnor) .lt. a) return z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) c*** We want y=UNI, then y1+s-2s*exp(-.5x^2), or ln(1+s-y)<-.5x^2+ln(2s), x=s*(b-x) if(alog(1.8857913-y) .lt. .5718733-.5*rnor**2) return 1 z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) v=4.656613e-10*(ishft(z,16)+iand(w,m2)) x=-alog(abs(v))/b z=36969*iand(z,m2)+ishft(z,-16) w=m1*iand(w,m2)+ishft(w,-16) y=.5+2.328306e-10*(ishft(z,16)+iand(w,m2)) y=-alog(y) if( y+y .lt. x*x) go to 1 rnor=sign(b+x,v) return end SUBROUTINE gaulag(x,w,n,alf) INTEGER n,MAXIT REAL alf,w(n),x(n) DOUBLE PRECISION EPS PARAMETER (EPS=3.D-14,MAXIT=10) CU USES gammln INTEGER i,its,j REAL ai,gammln DOUBLE PRECISION p1,p2,p3,pp,z,z1 do 13 i=1,n if(i.eq.1)then z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf) else if(i.eq.2)then z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n) else ai=i-2 z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))* *(z-x(i-2))/(1.+.3*alf) endif do 12 its=1,MAXIT p1=1.d0 p2=0.d0 do 11 j=1,n p3=p2 p2=p1 p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j 11 continue pp=(n*p1-(n+alf)*p2)/z z1=z z=z1-p1/pp if(abs(z-z1).le.EPS)goto 1 12 continue write(8,*)'Too many iterations in gaulag' 1 x(i)=z w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2) 13 continue return END