C:\win32app\salford\ftn77 spfit1 C:\win32app\salford\slink spfit1.obj -stack:0x20000000 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Driver for subroutine spfit c c Assume input files: c c file 1: raw data in n rows, 4 cols where c n: number of spatial data points c col 1: longitude c col 2: latitude c col 3: spatial variable c col 4: SD of spatial variable (0 if known exactly) c c file 2 (optional): correlations of W matrix c enter in order w(1,1), w(2,1), w(2,2), w(3,1), w(3,2), w(3,3), etc. c c file 3 (optional): additional covariates c n rows, nq columns, user prompted for nq c c file 4: coordinates for kriging output c matrix with npred rows, (nq+1) columns c col 1: longitude of kriged data point c col 2: latitude of kriged data point c cols 3 to nq+1: covariate of kriged data point c (NB it is assumed cthat first covariate is 1, i.e. intercept term) c c program also prompts for mod (form of covariance matrix) c and mdist (=0 for euclidean distance, =1 for geodesic) c parameter(n1max=100,n2max=5050,nqmax=10,nq2max=55,ndmax=3, 1 npredmax=873000,npmax=10) implicit double precision(a-h,o-z) dimension y(n1max),x(n1max,nqmax),s(n1max,ndmax),w(n2max), 1 v(n2max),b(npmax),bet(nqmax),vbet(nq2max),s0(npredmax,ndmax), 1 y0(npredmax),x0(npredmax,nqmax),vpred(npredmax), 1 w0(n1max),u1(n2max),x1(n1max,nqmax),y1(n1max),u2(nq2max), 1 h(npmax,npmax) dimension sd(n1max),corr(n2max) character*80 title common/lh/y,x,s,w,v,bet,vbet,u1,u2,x1,y1,n1,n2,nq,nq2,irem,mod, 1 nd,iw,ipro common /func1/ np common /dist/mdist write(*,*)' Enter name of main data matrix' read(*,'(a)')title open(7,file=title,status='old') nd=3 iw=1 n1=1 10 read(7,*,end=100)s(n1,1),s(n1,2),y(n1),sd(n1) s(n1,3)=1 n1=n1+1 if(n1.gt.n1max)then write(*,*)' Max observations reached' goto100 endif goto10 100 n1=n1-1 n2=n1*(n1+1)/2 close(unit=7) do 20 i=1,n2 corr(i)=0 20 continue do 30 i=1,n1 corr(i*(i+1)/2)=1 30 continue write(*,*)' Enter name of correlation matrix (0 if none)' read(*,'(a)')title if(title.eq.'0')goto200 open(7,file=title,status='old') do 40 i=1,n2 read(7,*)corr(i) 40 continue close(unit=7) 200 do 50 i=1,n1 do 50 j=1,i k=i*(i-1)/2+j w(k)=sd(i)*sd(j)*corr(k) 50 continue write(*,*)' Enter number of additional covariates' read(*,*)nq nq=nq+1 do 60 i=1,n1 x(i,1)=1 60 continue if(nq.eq.1)goto300 write(*,*)' Enter name of covariates file' read(*,'(a)')title open(7,file=title,status='old') do 70 i=1,n1 read(7,*)(x(i,j),j=2,nq) 70 continue close(unit=7) 300 continue write(*,*)' Enter kriging data file' read(*,'(a)')title open(7,file=title,status='old') npred=1 80 if(nq.eq.1)then read(7,*,end=400)s0(npred,1),s0(npred,2) else read(7,*,end=400)s0(npred,1),s0(npred,2),(x0(npred,j),j=2,nq) endif s0(npred,3)=1 x0(npred,1)=1 npred=npred+1 if(npred.gt.npredmax)then write(*,*)' Max kriging reached' goto400 endif goto80 400 continue close(unit=7) nq2=nq*(nq+1)/2 write(*,*)' Enter form of spatial covariance' write(*,*)' mod=0: nugget only (simple RE model)' write(*,*)' mod=1: Matern with nugget' write(*,*)' mod=2: exponential with nugget' c write(*,*)' mod=3: power law without nugget' c write(*,*)' mod=4: power law with nugget' write(*,*)' mod=5: Matern without nugget' write(*,*)' mod=6: exponential without nugget' write(*,*)' mod=7: Gaussian without nugget' read(*,*)mod if(mod.eq.0)np=1 if(mod.eq.1)np=4 if(mod.eq.2)np=3 if(mod.eq.3)np=3 if(mod.eq.4)np=4 if(mod.eq.5)np=3 if(mod.eq.6)np=2 if(mod.eq.7)np=2 write(*,*)' Enter distance type:' write(*,*)' 0 for Euclidean, 1 for geodesic' read(*,*)mdist c c test distance function c distmax=0 do 90 i=1,n1 do 90 j=1,i call dist1(s(i,1),s(i,2),s(j,1),s(j,2),h1) if(h1.gt.distmax)distmax=h1 90 continue write(*,*)' Largest distance is ',distmax write(*,*)' Largest log distance is ',dlog(distmax) write(*,*)' Enter starting parameters in order' if(mod.eq.0)write(*,*)' log nugget' if(mod.eq.1)write(*,*)' log scale, log range, kappa, log nugget' if(mod.eq.2)write(*,*)' log scale, log range, log nugget' if(mod.eq.5)write(*,*)' log scale, log range, kappa' if(mod.eq.6)write(*,*)' log scale, log range' if(mod.eq.7)write(*,*)' log scale, log range' read(*,*)(b(j),j=1,np) f0=func(b) write(*,*)' Initial value is ',f0 write(*,*)' Enter 1 to continue, 0 to stop' read(*,*)icont if(icont.eq.0)stop write(*,*)' Enter kriging output file' read(*,'(a)')title open(8,file=title,status='unknown') call spfit(b,f0,h,s0,x0,y0,vpred,npred,ih) write(*,1001)(b(j),j=1,np) write(*,1002)f0 do 110 i=1,npred-1 write(8,1001)s0(i,1),s0(i,2),y0(i),dsqrt(vpred(i)) 110 continue 1001 format(5f15.5) 1002 format(' NLLH = ',f40.15) end c c c subroutine spfit(y,x,s,w,v,b,bet,vbet,f0,s0,x0,y0,vpred,n1,n2, c 1 np,nq,nq2,nd,npred,irem) subroutine spfit(b,f0,h,s0,x0,y0,vpred,npred,ih) c c fit spatial model of form y = N[x'beta,v], then extend to kriging c c y: given vector of observations, dimension n1 bounded by n1max c x: given matrix of covariates, dimensions n1.nq, bounded by n1max.nqmax c s: given matrix of spatial locations, dimension n1.nd c ndmax: upper bound on nd c c w: given covariance matrix of measurement errors, of dimension c n2=n1*(n1+2)/2, expressed as a vector with entries w(1,1), w(1,2), w(2,2), w(3,1) c etc. n2 is bounded by n2max. If measurement error component is not required, c set w identical to 0. c c v: actual covariance matrix of y, expressed as a vector in same format as y c c b: unknown vector of parameters of the covariance matrix, of dimension np c at start, b has some (sensible) initial value; at end, b contains estimated c c h: unknown matrix of lead dimension ih c at end, h contains estimated inverse hessian matrix c c The user must supply a subroutine introduced by c subroutine cov(s1,s2,b,d0,cov,nd,np,mod,ifault) c which evaluates the covariance, as a function of given parameters b, between c two spatial locations defined by vectors s1 and s2. Here mod is used as an c indicator variable which could possibly index different parametric forms c (e.g. mod=1 means Matern, mod=2 means exponential, etc.) within a single "cov" c subroutine c c After evaluating v using subroutine cov, the program then adds the measurement c error component w to v c c On output bet contains the GLS estimators of the regression parameters beta, c and vbet is its covariance matrix (again written in vector format) c c irem=1 means REML estimate, irem=0 is MLE c c The function func(b) evaluates the neg log likelihood or REML equivalent c c f0: final value of func c c Kriging part: c c npred locations, maximum npredmax c c s0 defines locations of prediction points c x0 defines covariates at prediction points c y0 (output) are kriged values c vpred (output) are MSPEs c c Routines DFPMIN, LINMIN, MNBRAK, BRENT, F1DIM are copied from Numerical c Recipes with the following changes: c c 1. In all four routines insert the line c IMPLICIT DOUBLE PRECISION(A-H,O-Z) C after the opening line c c 2. Change the opening line of BRENT from c FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN) c to c DOUBLE PRECISION FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN) c c 3. Change the opening line of F1DIM from c PRECISION FUNCTION F1DIM(X) c to c DOUBLE PRECISION FUNCTION F1DIM(X) c c 4. Change all references to ABS( to DABS( c c 5. Change all references to SIGN( to DSIGN( c c 6. Change constants in PARAMETER statement: c c In DFPMIN change EPS=1.E-10 to EPS=1.D-14 c In LINMIN change TOL=1.E-4 to TOL=1.D-6 c In MNBRAK change TINY=1.E-20 to TINY=1.D-20 c In BRENT change ZEPS=1.0E-10 to ZEPS=1.0D-10 c c c these variables need to be defined as COMMON in the calling program: c c common/lh/y,x,s,w,v,bet,vbet,u1,u2,x1,y1,n1,n2,nq,nq2,irem,mod, c 1 nd,iw,ipro c common /func1/ np c c Note that various PARAMETER statements are used to bound the number of c variables: c n1max is max total number of observations c nqmax is max number of covariates c ndmax is max dimension of space c npredmax is max number of locations for prediction c npmax is max number of parameters in the covariance function c n2max must be n1max*(n1max+1)/2 c nq2max must be nqmax*(nqmax+1)/2 c c These can be changed, but both parameter(n1max=...) statements (as well c as any in the calling program) must be changed to the same values c parameter(n1max=100,n2max=5050,nqmax=10,nq2max=55,ndmax=3, 1 npredmax=873000,npmax=10) implicit double precision(a-h,o-z) dimension y(n1max),x(n1max,nqmax),s(n1max,ndmax),w(n2max), 1 v(n2max),b(npmax),bet(nqmax),vbet(nq2max),s0(npredmax,ndmax), 1 y0(npredmax),x0(npredmax,nqmax),vpred(npredmax), 1 w0(n1max),u1(n2max),x1(n1max,nqmax),y1(n1max),x00(nqmax), 1 u2(nq2max),s1(ndmax),s2(ndmax),w1(n1max),w2(nqmax),x2(nqmax), 1 h(npmax,npmax) common/lh/y,x,s,w,v,bet,vbet,u1,u2,x1,y1,n1,n2,nq,nq2,irem,mod, 1 nd,iw,ipro common /func1/ np common /wt/wt1,ztcov c c calculate MLE or REMLE c c if(ipro.eq.1.and.mod.eq.0)then c bb0=b(1) c b(1)=-9.99 c write(*,*)(b(j),j=1,np) c f0=func(b) c write(*,*)f0 c b(1)=bb0 c endif c write(*,*)np c write(*,*)(b(j),j=1,np) c f0=func(b) c write(*,*)f0 ftol=1.0d-6 iopt=1 c if(ipro.eq.1)then c write(*,*)' optimize? enter 1 for yes or no' c read(*,*)iopt c endif iopt=1 if(iopt.eq.1)then c call dfpmin(b,np,ftol,iter,f0) c write(*,*)' dfpmin done, f0= ',f0 call qnmin(np,b,h,f0,npmax,10000,ifail) write(*,*)' ifail = ',ifail endif c c set up for kriging c c c write(*,*)' check up on u1:' c do 201 i1=1,n1 c do 201 i4=1,n1 c vv=0 c do 202 i3=1,n1 c do 202 i2=1,n1 c if(i1.gt.i2)goto202 c if(i3.gt.i2)goto202 c k1=i2*(i2-1)/2+i1 c k2=i2*(i2-1)/2+i3 c k3=i3*(i3-1)/2+i4 c if(i4.gt.i3)k3=i4*(i4-1)/2+i3 c vv=vv+u1(k1)*u1(k2)*v(k3) c 202 continue c vv1=0 c if(i1.eq.i4)vv1=1 c if(dabs(vv-vv1).gt.1.0d-10)write(*,*)i1,i2,vv c 201 continue c write(*,*)' press return to continue' c read(*,'(a)')title c return c do 10 ipred=1,npred do 20 j=1,nd s1(j)=s0(ipred,j) s2(j)=s0(ipred,j) 20 continue call cov(s1,s2,b,d0,v0,nd,np,mod,ifault) do 40 i=1,n1 do 30 j=1,nd s1(j)=s(i,j) s2(j)=s0(ipred,j) 30 continue call cov(s1,s2,b,d0,w0(i),nd,np,mod,ifault) 40 continue do 50 j=1,nq x00(j)=x0(ipred,j) 50 continue call krig(y,x,v,y0(ipred),x00,vpred(ipred),v0,w0,bet,u1,u2,x1,y1, 1 w1,w2,x2,n1max,n2max,nqmax,nq2max,n1,n2,nq,nq2) if(20000*(ipred/20000).eq.ipred)write(*,*)' ipred = ',ipred 10 continue return end c subroutine krig(y,x,v,y0,x0,vpred,v0,w0,bet,u1,u2,x1,y1, 1 w1,w2,x2,n1max,n2max,nqmax,nq2max,n1,n2,nq,nq2) implicit double precision(a-h,o-z) dimension y(n1max),x(n1max,nqmax),v(n2max),x0(nqmax),w0(n1max), 1 bet(nqmax),u1(n2max),x1(n1max,nqmax),y1(n1max),u2(nq2max), 1 w1(n1max),w2(nqmax),x2(nqmax) c c Kriging routine c c y: observed vector, dimension n1 c x: covariate matrix of observations, dimension n1.nq c v: covariance matrix of observations, in vector format c c y0: desired predictions (on output contains kriged value) c x0: covariate vector associated with predictions, dimension nq c w0: cross-covariance matrix, dimension n1 c v0: unconditional variance of y0 (input) c vpred: prediction variance of y0 (output) c c bet, u1, u2, x1, y1 are input and assumed to have been calculated by c a previous evaluation of func c c Step 1: Compute w1=u1'w0 c do 10 i1=1,n1 w1(i1)=0 do 10 i2=1,i1 k=i1*(i1-1)/2+i2 w1(i1)=w1(i1)+u1(k)*w0(i2) 10 continue c c Step 2: Compute sum of squares of w_1 (same as $w_0^TV^{-1}w_0$) c - starts computation of vpred c vpred=v0 do 20 i=1,n1 vpred=vpred-w1(i)**2 20 continue c c Step 3: Compute x2=x0-x1'w1 c Step 4: Compute w2=u2'x2 c do 30 j=1,nq x2(j)=x0(j) do 30 i=1,n1 x2(j)=x2(j)-x1(i,j)*w1(i) 30 continue do 40 j1=1,nq w2(j1)=0 do 40 j2=1,j1 k=j1*(j1-1)/2+j2 w2(j1)=w2(j1)+u2(k)*x2(j2) 40 continue c c Step 5: Compute y0=w1'y1+x2'bet c y0=0 do 60 i=1,n1 y0=y0+w1(i)*y1(i) 60 continue do 70 j=1,nq y0=y0+x2(j)*bet(j) 70 continue c c Complete computation of vpred c do 80 j=1,nq vpred=vpred+w2(j)*w2(j) 80 continue return end c double precision function func(b) implicit double precision(a-h,o-z) parameter(n1max=100,n2max=5050,nqmax=10,nq2max=55,ndmax=3, 1 npredmax=873000,npmax=10) dimension y(n1max),x(n1max,nqmax),s(n1max,ndmax),w(n2max), 1 v(n2max),b(npmax),bet(nqmax),vbet(nq2max),u1(n2max), 1 x1(n1max,nqmax),y1(n1max),y2(nqmax),x2(nq2max),u2(nq2max), 1 y3(nqmax),s1(ndmax),s2(ndmax) common/lh/y,x,s,w,v,bet,vbet,u1,u2,x1,y1,n1,n2,nq,nq2,irem,mod, 1 nd,iw,ipro common /func1/ np data pi/3.141592654d0/ c c Evaluate covariance function and incorporate w c c use ifault as an indicator that parameters are infeasible c c write(*,1001)(b(j),j=1,np) c if(b(np).lt.1.or.b(np).gt.10)goto1000 c write(*,*)' n1 = ',n1 do 10 i1=1,n1 do 10 i2=1,i1 c write(*,*)' i1, i2',i1,i2 k=i1*(i1-1)/2+i2 do 15 j=1,nd s1(j)=s(i1,j) s2(j)=s(i2,j) 15 continue call cov(s1,s2,b,d0,cov1,nd,np,mod,ifault) if(ifault.eq.1)goto1000 v(k)=cov1 if(iw.eq.1.or.(iw.eq.2.and.i1.eq.i2))v(k)=cov1+w(k) c distmax=exp(b(np)) c if(d0.gt.distmax)then c v(k)=cov1 c else c d5=d0/distmax c v(k)=cov1+2*w(k)*(dacos(d5)-d5*dsqrt(1-d5*d5))/pi c endif 10 continue c c Step 1: write v=u1'u1 where u1 is upper trangular c call chol(v,n1,n2,u1,nullty,ifault) if(nullty.ne.0)write(*,*)' chol: nullty = ',nullty if(ifault.ne.0)write(*,*)' chol: ifault = ',ifault if(ifault.ne.0.or.nullty.ne.0)goto1000 c c Step 2: Replace u1 by its inverse (also upper triangular) c call cholinv(u1,n1,n2) c c Step 3: Calculate x1=u1'x c do 30 i1=1,n1 do 30 j=1,nq x1(i1,j)=0 do 30 i2=1,i1 k=i1*(i1-1)/2+i2 x1(i1,j)=x1(i1,j)+u1(k)*x(i2,j) 30 continue c c Step 4: Calculate y1=u1'y. Also y1sq is sum of squares of y1 c do 40 i1=1,n1 y1(i1)=0 do 40 i2=1,i1 k=i1*(i1-1)/2+i2 y1(i1)=y1(i1)+u1(k)*y(i2) 40 continue y1sq=0 do 50 i=1,n1 y1sq=y1sq+y1(i)*y1(i) 50 continue c write(*,*)' step @4' c c Step 5: Compute y2=x1'y1 (this is $X^TV^{-1}Y$) c do 60 j=1,nq y2(j)=0 do 60 i=1,n1 y2(j)=y2(j)+x1(i,j)*y1(i) 60 continue c write(*,*)' step @5' c c Step 6: Compute x2=x1'x1 in vector form (this is $X^TV^{-1}X$) c do 70 j1=1,nq do 70 j2=1,j1 k=j1*(j1-1)/2+j2 x2(k)=0 do 70 i=1,n1 x2(k)=x2(k)+x1(i,j1)*x1(i,j2) 70 continue c write(*,*)' step @6' c c Step 7: u2 is Cholesky decomp of x2 - then invert u2 c call chol(x2,nq,nq2,u2,nullty,ifault) if(nullty.ne.0)write(*,*)' chol: nullty = ',nullty if(ifault.ne.0)write(*,*)' chol: ifault = ',ifault if(ifault.ne.0.or.nullty.ne.0)goto1000 call cholinv(u2,nq,nq2) c write(*,*)' step @7' c c Step 8: Sum log diagonal entries of u1 - this is $-{1\over 2}\log |V|$ c Step 9: Sum log diagonal entries of u2 - this is $-{1\over 2}\log |X^TV^{-1}X|$ c detlog1=0 detlog2=0 do 80 i=1,n1 k=i*(i+1)/2 detlog1=detlog1+dlog(u1(k)) 80 continue do 90 i=1,nq k=i*(i+1)/2 detlog2=detlog2+dlog(u2(k)) 90 continue c write(*,*)' step @9' c c Step 10: calculate vbet=u2.u2' c do 100 j1=1,nq do 100 j2=1,j1 k1=j1*(j1-1)/2+j2 vbet(k1)=0 do 100 j3=j1,nq k2=j3*(j3-1)/2+j1 k3=j3*(j3-1)/2+j2 vbet(k1)=vbet(k1)+u2(k2)*u2(k3) 100 continue c write(*,*)' step @10' c c Step 11: Calculate y3=u2'y2 and sum of squares y3sq c do 110 j1=1,nq y3(j1)=0 do 110 j2=1,j1 k=j1*(j1-1)/2+j2 y3(j1)=y3(j1)+u2(k)*y2(j2) 110 continue y3sq=0 do 120 j=1,nq y3sq=y3sq+y3(j)*y3(j) 120 continue c write(*,*)' step @11' c c Step 12: Calculate bet=u2.y3 c do 130 j1=1,nq bet(j1)=0 do 130 j2=j1,nq k=j2*(j2-1)/2+j1 bet(j1)=bet(j1)+u2(k)*y3(j2) 130 continue c write(*,*)' step @12' c c Step 13: Calculate negative (restricted) log likelihood c func=0.5*(y1sq-y3sq)-detlog1-irem*detlog2 c write(*,1001)(b(j),j=1,np) c write(*,1002)func c write(*,1003)y1sq-y3sq,detlog1,detlog2 return 1000 func=1.0d10 c write(*,1001)(b(j),j=1,np) c write(*,1002)func 1001 format(10f10.5) 1002 format(' NLLH = ',f35.10) c 1003 format(' G2, detlog1, detlog2 ',3f15.8) return end c subroutine dfunc(b,g1) implicit double precision(a-h,o-z) parameter(npmax=10) dimension b(npmax),g1(npmax) common /func1/ np data eps/1.0d-6/ f0=func(b) do 10 i=1,np b(i)=b(i)+eps f1=func(b) b(i)=b(i)-eps g1(i)=(f1-f0)/eps 10 continue return end c c c subroutine cov(s1,s2,b,d0,cov1,nd,np,mod,ifault) c c Compute covariance two spatial locations defined by vectors s1 and s2 c c This is a dummy routine to allow calculation of spatial averages to be c handled by the program c c If s1(3) and s2(3) are both 1, treat as a 2-D vector and go to cov2 c Otherwise treat as average c parameter(n1max=100,ndmax=3) implicit double precision(a-h,o-z) dimension s1(nd),s2(nd),vv(1),b(np),dist(1),theta(3) dimension scov(n1max,ndmax),ipop(200) common /icov/ scov,ipop,ptot,ncov if(s1(3).lt.1.5.and.s2(3).lt.1.5)then call cov2(s1,s2,b,d0,cov1,nd,np,mod,ifault) return endif if(s1(3).gt.1.5d0.and.s2(3).lt.1.5d0)goto100 if(s1(3).lt.1.5d0.and.s2(3).gt.1.5d0)goto200 sum1=0 do 10 i=1,ncov do 10 j=1,ncov s1(1)=scov(i,1) s1(2)=scov(i,2) s2(1)=scov(j,1) s2(2)=scov(j,2) call cov2(s1,s2,b,d0,cov1,nd,np,mod,ifault) if(ifault.eq.1)return sum1=sum1+cov1*ipop(i)*ipop(j) 10 continue cov1=sum1/ptot**2 c write(*,*)' cov1 = ',cov1 return 100 sum1=0 do 110 i=1,ncov s1(1)=scov(i,1) s1(2)=scov(i,2) call cov2(s1,s2,b,d0,cov1,nd,np,mod,ifault) if(ifault.eq.1)return sum1=sum1+cov1*ipop(i) 110 continue cov1=sum1/ptot c write(*,*)' Covariance ',cov1 return 200 sum1=0 do 210 i=1,ncov s2(1)=scov(i,1) s2(2)=scov(i,2) call cov2(s1,s2,b,d0,cov1,nd,np,mod,ifault) if(ifault.eq.1)return sum1=sum1+cov1*ipop(i) 210 continue cov1=sum1/ptot c write(*,*)' Covariance2 ',cov1 return end c c c subroutine cov2(s1,s2,b,d0,cov1,nd,np,mod,ifault) c c Actual covariance calculation c c mod is an indicator allowing different models to be evaluated c within the same routine c c nd is dimension of space c c this version assumes matern covariances with geodesic distances c assumes first two coordinates represent latitude and longitude respectively c c np is number of parameters, represented by vector b c c d0 is distance between s1 and s2 c c ifault is an error indicator - use to catch inappropriate parameter values c implicit double precision(a-h,o-z) dimension s1(nd),s2(nd),vv(1),b(np),dist(1),theta(3) ifault=0 call dist1(s1(1),s1(2),s2(1),s2(2),dist(1)) d0=dist(1) if(mod.eq.0)goto60 if(mod.eq.2)goto20 if(mod.eq.5)goto50 if(mod.eq.6)goto70 if(mod.eq.7)goto80 if(mod.gt.2)goto30 c c mod=1 - Matern with nugget c if(dabs(b(1)).gt.10)ifault=1 if(dabs(b(2)).gt.10)ifault=1 if(dabs(b(4)).gt.10)ifault=1 if(b(2).lt.0)ifault=1 if(b(3).lt.0.001)ifault=1 if(ifault.eq.1)return theta(1)=exp(b(1)) theta(2)=exp(b(2)) theta(3)=b(3) if(b(3).gt.50)theta(3)=50 call rkmat(theta,dist,vv,1) cov1=vv(1)+exp(b(4)) return c c mod=2 - exponential with nugget c 20 if(dabs(b(1)).gt.20)ifault=1 if(dabs(b(2)).gt.20)ifault=1 if(dabs(b(3)).gt.20)ifault=1 if(b(2).lt.0)ifault=1 if(ifault.eq.1)return theta(1)=exp(b(1)) theta(2)=exp(b(2)) cov1=theta(1)*exp(-dist(1)/theta(2)) if(dist(1).lt.1.0d-10)cov1=cov1+exp(b(3)) return c c mod=3 or 4 - power law model (mod=4 means with nugget) c 30 if(dabs(b(1)).gt.20)ifault=1 if(b(2).gt.0.693147.or.b(2).lt.-7)ifault=1 if(mod.eq.4.and.dabs(b(3)).gt.20)ifault=1 if(ifault.eq.1)return b0=0 if(mod.eq.4)b0=exp(b(3)) cov1=b0+exp(b(1))*15000**exp(b(2)) if(d0.gt.0.001)cov1=cov1-b0-exp(b(1))*d0**exp(b(2)) return c c mod=5 - Matern without nugget c 50 if(dabs(b(1)).gt.10)ifault=1 if(dabs(b(2)).gt.10)ifault=1 if(b(2).lt.0)ifault=1 if(b(3).lt.0.001)ifault=1 if(ifault.eq.1)return theta(1)=exp(b(1)) theta(2)=exp(b(2)) theta(3)=b(3) if(b(3).gt.50)theta(3)=50 call rkmat(theta,dist,vv,1) cov1=vv(1) return c c mod=0 nugget only c 60 if(dabs(b(1)).gt.8)then ifault=1 return endif if(dist(1).gt.1.0d-5)then cov1=0 else cov1=exp(b(1)) endif return c c mod=6 exponential without nugget c 70 if(dabs(b(1)).gt.20)ifault=1 if(dabs(b(2)).gt.20)ifault=1 if(b(2).lt.0)ifault=1 if(ifault.eq.1)return theta(1)=exp(b(1)) theta(2)=exp(b(2)) cov1=theta(1)*exp(-dist(1)/theta(2)) return c c mod=7 gaussian without nugget c 80 if(dabs(b(1)).gt.20)ifault=1 if(dabs(b(2)).gt.20)ifault=1 if(b(2).lt.0)ifault=1 if(ifault.eq.1)return theta(1)=exp(b(1)) theta(2)=exp(b(2)) cov1=theta(1)*exp(-(dist(1)/theta(2))**2) return end c SUBROUTINE DIST1(X1,Y1,X2,Y2,H) IMPLICIT DOUBLE PRECISION(A-H,O-Z) common /dist/mdist DATA PI/3.141592654D0/ if(mdist.eq.0)then d2=dsqrt((x1-x2)**2+(y1-y2)**2) return endif c c calculate geodesic distances in km: c x1,x2 are longitudes, y1, y2 are latitudes c Modify to allow for geodesic distances, in units of 100km D2=(DCOS(Y1*PI/180)*DCOS(X1*PI/180) 1 -DCOS(Y2*PI/180)*DCOS(X2*PI/180))**2 1 +(DCOS(Y1*PI/180)*DSIN(X1*PI/180) 1 -DCOS(Y2*PI/180)*DSIN(X2*PI/180))**2 1 +(DSIN(Y1*PI/180)-DSIN(Y2*PI/180))**2 H=12732.4*DASIN(DSQRT(D2)/2) RETURN END c C SUBROUTINE QNMIN(N,B,H,P0,IH,NEVALS,IFAIL) C---------------------------------------------------------- C QUASI-NEWTON MINIMIZER FROM NASH, "COMPACT NUMERICAL C METHODS FOR COMPUTERS" (1979), ALG 21. C C B.D. RIPLEY 5/1980 C C ADAPTED BY R.L. SMITH 4/1988 C C CALLS ROUTINES C FUNCT1(N,B,P0) : RETURNS VALUE P0 AT C PARAMETER VECTOR B OF DIMENSION N, C GRAD(N,B,G,P0) : FOLLOWS FUNCT1 (P0 INPUT) C AND GIVES GRADIENT VECTOR IN G. C C H IS FINAL ESTIMATED INVERSE HESSIAN MATRIX OF LEAD C DIMENSION IH. C C NEVALS: TOTAL NUMBER OF FUNCTION EVALS. PERMITTED C (SET BY USER) C IFAIL: FAILURE INDICATOR (=0 MEANS OK) C (PART OF OUTPUT FROM PROGRAM) C C IF THE INITIAL FUNCTION EVALUATION IS LARGER THAN C 1.0D9 THE PROGRAM EXITS AT ONCE - THIS SHOULD BE C USED AS AN INDICATOR OF AN INFEASIBLE POINT. C---------------------------------------------------------- IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION B(N), H(IH,N) DIMENSION X(23), C(23), G(23), T(23), IBCD(7) DOUBLE PRECISION K INTEGER COUNT DATA W,TOL/0.2,1.0D-4/,EPS/1.0D-6/ IF (N.LT.0.OR.N.GT.23) GO TO 160 IFN = N+1 IG = 1 P0=FUNC(B) C CALL FUNCT1(N,B,P0) IF(P0.GT.1.0D9)GOTO180 CALL DFUNC(B,G) C CALL GRAD(N,B,G,P0) C C RESET HESSIAN C 10 DO 30 I = 1,N DO 20 J = 1,N 20 H(I,J) = 0.0 30 H(I,I) = 1.0 ILAST = IG C C TOP OF ITERATION C 40 DO 50 I = 1,N X(I) = B(I) 50 C(I) = G(I) C C FIND SEARCH DIRECTION T C D1 = 0.0 SN=0.0 DO 70 I = 1,N S = 0.0 DO 60 J = 1,N 60 S = S-H(I,J)*G(J) T(I) = S SN = SN+S*S 70 D1 = D1-S*G(I) C C CHECK IF DOWNHILL C IF (D1.LE.0.0.AND.ILAST.NE.IG) GO TO 10 IF (D1.LE.0.0.AND.ILAST.EQ.IG) GO TO 150 C C SEARCH ALONG T C SN = 0.5/DSQRT(SN) K = DMIN1(1.0D0,SN) 80 COUNT = 0 DO 90 I = 1,N B(I) = X(I)+K*T(I) IF (DABS(B(I)-X(I)).LT.EPS) COUNT = COUNT+1 90 CONTINUE C C CHECK IF CONVERGED C IF (COUNT.EQ.N) GO TO 150 C CALL FUNCT1(N,B,P) P=FUNC(B) IFN = IFN+1 IF (IFN.GE.NEVALS) GO TO 170 IF (P.LT.P0-D1*K*TOL) GO TO 100 K = W*K GO TO 80 C C NEW LOWEST VALUE C 100 P0 = P IG = IG+1 C CALL GRAD(N,B,G,P) CALL DFUNC(B,G) IFN = IFN+N C C UPDATE HESSIAN C D1 = 0.0 DO 110 I = 1,N T(I) = K*T(I) C(I) = G(I)-C(I) 110 D1 = D1+T(I)*C(I) C C CHECK IF +VE DEF ADDITION C IF (D1.LE.0.0) GO TO 10 D2 = 0.0 DO 130 I = 1,N S = 0.0 DO 120 J = 1,N 120 S = S+H(I,J)*C(J) X(I) = S 130 D2 = D2+S*C(I) D2 = 1+D2/D1 DO 140 I = 1,N DO 140 J = 1,N 140 H(I,J) = H(I,J)-(T(I)*X(J)+T(J)*X(I)-D2*T(I)*T(J))/D1 GO TO 40 150 IFAIL = 0 C SUCCESSFUL CONCLUSION RETURN 160 IFAIL = 1 C N OUT OF RANGE RETURN 170 IFAIL = 2 C TOO MANY FUNCTION EVALUATIONS RETURN 180 IFAIL = 3 C INITIAL POINT INFEASIBLE RETURN END C c SUBROUTINE CHOLINV(C,N,NN) C C FOLLOWS CHOL - INVERT LOWER TRIANGULAR MATRIX CONTAINED IN C C IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION C(NN) C(1)=1/C(1) DO 10 I=2,N DO 20 J=1,I-1 SUM=0 DO 30 K=J,I-1 30 SUM=SUM+C(I*(I-1)/2+K)*C(K*(K-1)/2+J) 20 C(I*(I-1)/2+J)=-SUM/C(I*(I+1)/2) 10 C(I*(I+1)/2)=1/C(I*(I+1)/2) RETURN END C SUBROUTINE CHOL(A,N,NN,U,NULLTY,IFAULT) C C ALGORITHM AS 6 APPL. STATIST. (1968) VOL. 17, P.195 C C GIVEN A SYMMETRIC MATRIC ORDER N AS LOWER TRIANGLE IN A( ) C CALCULATES AN UPPER TRIANGLE, U( ), SUCH THAT UPRIME * U = A. C A MUST BE POSITIVE SEMI-DEFINITE. ETA IS SET TO MULTIPLYING C FACTOR DETERMINING EFFECTIVE ZERO FOR PIVOT C C A: ARRAY(NN), THE INPUT MATRIX, STORED IN ORDER A(1,1), C A(2,1), A(2,2), A(3,1), A(3,2), A(3,3),... C C N: INTEGER, INPUT, THE ORDER OF A C C NN: INTEGER, INPUT, MUST BE N*(N+1)/2 C C U: ARRAY(NN), OUTPUT, IN ORDER U(1,1), U(1,2), U(2,2), C U(1,3), U(2,3), U(3,3),... C C NULLTY: INTEGER, OUTPUT, THE NULLITY OF A C C IFAULT: INTEGER, OUTPUT, =1 IF N IS <1, =2 IF A NOT POSITIVE SD, C =3 IF NN.NE.N*(N+1)/2, 0 OTHERWISE C DOUBLE PRECISION A(NN),U(NN),ETA,ETA2,X,W,ZERO,ZABS,ZSQRT C C THE VALUE OF ETA WILL DEPEND ON THE WORD-LENGTH OF THE COMPUTER C BEING USED. SUGGEST 1E-5 FOR SINGLE PRECISION, 1D-9 FOR D.P. C DATA ETA, ZERO/1.0D-9, 0.0/ C ZABS(X)=DABS(X) ZSQRT(X)=DSQRT(X) C IFAULT=1 IF(N.LE.0)RETURN IFAULT=3 IF(NN.NE.N*(N+1)/2)RETURN IFAULT=2 NULLTY=0 J=1 K=0 ETA2=ETA*ETA II=0 DO 80 ICOL=1,N II=II+ICOL X=ETA2*A(II) L=0 KK=0 DO 40 IROW=1,ICOL KK=KK+IROW K=K+1 W=A(K) M=J DO 10 I=1,IROW L=L+1 IF(I.EQ.IROW)GOTO20 W=W-U(L)*U(M) M=M+1 10 CONTINUE 20 IF(IROW.EQ.ICOL)GOTO50 IF(U(L).EQ.ZERO)GOTO30 U(K)=W/U(L) GOTO40 30 IF(W*W.GT.ZABS(ETA*A(KK)))RETURN U(K)=ZERO 40 CONTINUE 50 IF(ZABS(W).LE.ZABS(ETA*A(K)))GOTO60 IF(W.LT.ZERO)RETURN U(K)=ZSQRT(W) GOTO70 60 U(K)=ZERO NULLTY=NULLTY+1 70 J=J+ICOL 80 CONTINUE IFAULT=0 RETURN END C c subroutine rkmat(theta,x,y,n) c c msh: date: 1/13/89 c c purpose: calculate Matern class covariances c If the smoothness is above 50 the Squared c Exponential limiting covariance is used. c c arguments: c theta parameter vector (scale, range, smoothness) c x vector of distances (input) c y vector of Matern covariances (output) c n length of x (and y) c c for a technical introduction see c c Handcock and Stein (1993), c 'A Bayesian Analysis of Kriging', c Technometrics, 35, 4, 403-410. c c or contact c c mhandcock@stern.nyu.edu c c amendments: c c date description c ---- ----------- c c implicit double precision (a-h,o-z) parameter( zero = 0.0d0, one = 1.0d0 ) parameter( half = 0.5d0, two = 2.0d0 ) parameter( fifty = 50.0d0 ) c dimension theta(3),x(n),y(n) dimension bk(50) c t1 = theta(1) t2l = one/theta(2) t3 = theta(3) c t2 = two*dsqrt(t3)*t2l c do 10 i = 1, n d = x(i) if( d .le. zero ) then y(i) = t1 else if(t3 .lt. fifty) then cnv = (two**(t3-one))*dgamma(t3) p1 = t2 * d it3 = dint(t3) call rkbesl(p1,t3-it3,it3+1,1,bk,ncalc) y(i) = t1*(p1**t3)*bk(it3+1)/cnv else p1 = t2l * d y(i) = t1*dexp(-p1*p1) endif endif 10 continue c c finish up c return end c SUBROUTINE RKBESL(X,ALPHA,NB,IZE,BK,NCALC) C------------------------------------------------------------------- C C This FORTRAN 77 routine calculates modified Bessel functions C of the second kind, K SUB(N+ALPHA) (X), for non-negative C argument X, and non-negative order N+ALPHA, with or without C exponential scaling. C C Explanation of variables in the calling sequence C C Description of output values .. C C X - Working precision non-negative real argument for which C K's or exponentially scaled K's (K*EXP(X)) C are to be calculated. If K's are to be calculated, C X must not be greater than XMAX (see below). C ALPHA - Working precision fractional part of order for which C K's or exponentially scaled K's (K*EXP(X)) are C to be calculated. 0 .LE. ALPHA .LT. 1.0. C NB - Integer number of functions to be calculated, NB .GT. 0. C The first function calculated is of order ALPHA, and the C last is of order (NB - 1 + ALPHA). C IZE - Integer type. IZE = 1 if unscaled K's are to be calculated, C and 2 if exponentially scaled K's are to be calculated. C BK - Working precision output vector of length NB. If the C routine terminates normally (NCALC=NB), the vector BK C contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), C or the corresponding exponentially scaled functions. C If (0 .LT. NCALC .LT. NB), BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C NCALC - Integer output variable indicating possible errors. C Before using the vector BK, the user should check that C NCALC=NB, i.e., all orders have been calculated to C the desired accuracy. See error returns below. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta = Radix for the floating-point system C minexp = Smallest representable power of beta C maxexp = Smallest power of beta that overflows C EPS = The smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMAX = Upper limit on the magnitude of X when IZE=1; Solution C to equation: C W(X) * (1-1/8X+9/128X**2) = beta**minexp C where W(X) = EXP(-X)*SQRT(PI/2X) C SQXMIN = Square root of beta**minexp C XINF = Largest positive machine number; approximately C beta**maxexp C XMIN = Smallest positive machine number; approximately C beta**minexp C C C Approximate values for some important machines are: C C beta minexp maxexp EPS C C CRAY-1 (S.P.) 2 -8193 8191 7.11E-15 C Cyber 180/185 C under NOS (S.P.) 2 -975 1070 3.55E-15 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 -126 128 1.19E-7 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 -1022 1024 2.22D-16 C IBM 3033 (D.P.) 16 -65 63 2.22D-16 C VAX (S.P.) 2 -128 127 5.96E-8 C VAX D-Format (D.P.) 2 -128 127 1.39D-17 C VAX G-Format (D.P.) 2 -1024 1023 1.11D-16 C C C SQXMIN XINF XMIN XMAX C C CRAY-1 (S.P.) 6.77E-1234 5.45E+2465 4.59E-2467 5674.858 C Cyber 180/855 C under NOS (S.P.) 1.77E-147 1.26E+322 3.14E-294 672.788 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 1.08E-19 3.40E+38 1.18E-38 85.337 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.49D-154 1.79D+308 2.23D-308 705.342 C IBM 3033 (D.P.) 7.35D-40 7.23D+75 5.40D-79 177.852 C VAX (S.P.) 5.42E-20 1.70E+38 2.94E-39 86.715 C VAX D-Format (D.P.) 5.42D-20 1.70D+38 2.94D-39 86.715 C VAX G-Format (D.P.) 7.46D-155 8.98D+307 5.57D-309 706.728 C C******************************************************************* C******************************************************************* C C Error returns C C In case of an error, NCALC .NE. NB, and not all K's are C calculated to the desired accuracy. C C NCALC .LT. -1: An argument is out of range. For example, C NB .LE. 0, IZE is not 1 or 2, or IZE=1 and ABS(X) .GE. C XMAX. In this case, the B-vector is not calculated, C and NCALC is set to MIN0(NB,0)-2 so that NCALC .NE. NB. C NCALC = -1: Either K(ALPHA,X) .GE. XINF or C K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) .GE. XINF. In this case, C the B-vector is not calculated. Note that again C NCALC .NE. NB. C C 0 .LT. NCALC .LT. NB: Not all requested function values could C be calculated accurately. BK(I) contains correct function C values for I .LE. NCALC, and contains the ratios C K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. C C C Intrinsic functions required are: C C ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT C C C Acknowledgement C C This program is based on a program written by J. B. Campbell C (2) that computes values of the Bessel functions K of real C argument and real order. Modifications include the addition C of non-scaled functions, parameterization of machine C dependencies, and the use of more accurate approximations C for SINH and SIN. C C References: "On Temme's Algorithm for the Modified Bessel C Functions of the Third Kind," Campbell, J. B., C TOMS 6(4), Dec. 1980, pp. 581-586. C C "A FORTRAN IV Subroutine for the Modified Bessel C Functions of the Third Kind of Real Order and Real C Argument," Campbell, J. B., Report NRC/ERB-925, C National Research Council, Canada. C C Latest modification: May 30, 1989 C C Modified by: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C------------------------------------------------------------------- INTEGER I,IEND,ITEMP,IZE,J,K,M,MPLUS1,NB,NCALC DOUBLE PRECISION 1 A,ALPHA,BLPHA,BK,BK1,BK2,C,D,DM,D1,D2,D3,ENU,EPS,ESTF,ESTM, 2 EX,FOUR,F0,F1,F2,HALF,ONE,P,P0,Q,Q0,R,RATIO,S,SQXMIN,T,TINYX, 3 TWO,TWONU,TWOX,T1,T2,WMINF,X,XINF,XMAX,XMIN,X2BY4,ZERO DIMENSION BK(1),P(8),Q(7),R(5),S(4),T(6),ESTM(6),ESTF(7) C--------------------------------------------------------------------- C Mathematical constants C A = LOG(2.D0) - Euler's constant C D = SQRT(2.D0/PI) C--------------------------------------------------------------------- DATA HALF,ONE,TWO,ZERO/0.5D0,1.0D0,2.0D0,0.0D0/ DATA FOUR,TINYX/4.0D0,1.0D-10/ DATA A/ 0.11593151565841244881D0/,D/0.797884560802865364D0/ C--------------------------------------------------------------------- C Machine dependent parameters C--------------------------------------------------------------------- DATA EPS/2.22D-16/,SQXMIN/1.49D-154/,XINF/1.79D+308/ DATA XMIN/2.23D-308/,XMAX/705.342D0/ C--------------------------------------------------------------------- C P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA C + Euler's constant C Coefficients converted from hex to decimal and modified C by W. J. Cody, 2/26/82 C R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) C T - Approximation for SINH(Y)/Y C--------------------------------------------------------------------- DATA P/ 0.805629875690432845D00, 0.204045500205365151D02, 1 0.157705605106676174D03, 0.536671116469207504D03, 2 0.900382759291288778D03, 0.730923886650660393D03, 3 0.229299301509425145D03, 0.822467033424113231D00/ DATA Q/ 0.294601986247850434D02, 0.277577868510221208D03, 1 0.120670325591027438D04, 0.276291444159791519D04, 2 0.344374050506564618D04, 0.221063190113378647D04, 3 0.572267338359892221D03/ DATA R/-0.48672575865218401848D+0, 0.13079485869097804016D+2, 1 -0.10196490580880537526D+3, 0.34765409106507813131D+3, 2 0.34958981245219347820D-3/ DATA S/-0.25579105509976461286D+2, 0.21257260432226544008D+3, 1 -0.61069018684944109624D+3, 0.42269668805777760407D+3/ DATA T/ 0.16125990452916363814D-9, 0.25051878502858255354D-7, 1 0.27557319615147964774D-5, 0.19841269840928373686D-3, 2 0.83333333333334751799D-2, 0.16666666666666666446D+0/ DATA ESTM/5.20583D1, 5.7607D0, 2.7782D0, 1.44303D1, 1.853004D2, 1 9.3715D0/ DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0, 1 8.45096D1, 2.0D1/ C--------------------------------------------------------------------- EX = X ENU = ALPHA NCALC = MIN(NB,0)-2 IF ((NB .GT. 0) .AND. ((ENU .GE. ZERO) .AND. (ENU .LT. ONE)) 1 .AND. ((IZE .GE. 1) .AND. (IZE .LE. 2)) .AND. 2 ((IZE .NE. 1) .OR. (EX .LE. XMAX)) .AND. 3 (EX .GT. ZERO)) THEN K = 0 IF (ENU .LT. SQXMIN) ENU = ZERO IF (ENU .GT. HALF) THEN K = 1 ENU = ENU - ONE END IF TWONU = ENU+ENU IEND = NB+K-1 C = ENU*ENU D3 = -C IF (EX .LE. ONE) THEN C--------------------------------------------------------------------- C Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA C Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA C--------------------------------------------------------------------- D1 = ZERO D2 = P(1) T1 = ONE T2 = Q(1) DO 10 I = 2,7,2 D1 = C*D1+P(I) D2 = C*D2+P(I+1) T1 = C*T1+Q(I) T2 = C*T2+Q(I+1) 10 CONTINUE D1 = ENU*D1 T1 = ENU*T1 F1 = LOG(EX) F0 = A+ENU*(P(8)-ENU*(D1+D2)/(T1+T2))-F1 Q0 = EXP(-ENU*(A-ENU*(P(8)+ENU*(D1-D2)/(T1-T2))-F1)) F1 = ENU*F0 P0 = EXP(F1) C--------------------------------------------------------------------- C Calculation of F0 = C--------------------------------------------------------------------- D1 = R(5) T1 = ONE DO 20 I = 1,4 D1 = C*D1+R(I) T1 = C*T1+S(I) 20 CONTINUE IF (ABS(F1) .LE. HALF) THEN F1 = F1*F1 D2 = ZERO DO 30 I = 1,6 D2 = F1*D2+T(I) 30 CONTINUE D2 = F0+F0*F1*D2 ELSE D2 = SINH(F1)/ENU END IF F0 = D2-ENU*D1/(T1*P0) IF (EX .LE. TINYX) THEN C-------------------------------------------------------------------- C X.LE.1.0E-10 C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK(1) = F0+EX*F0 IF (IZE .EQ. 1) BK(1) = BK(1)-EX*BK(1) RATIO = P0/F0 C = EX*XINF IF (K .NE. 0) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X), C ALPHA .GE. 1/2 C-------------------------------------------------------------------- NCALC = -1 IF (BK(1) .GE. C/RATIO) GO TO 500 BK(1) = RATIO*BK(1)/EX TWONU = TWONU+TWO RATIO = TWONU END IF NCALC = 1 IF (NB .EQ. 1) GO TO 500 C-------------------------------------------------------------------- C Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), L = 1, 2, ... , NB-1 C-------------------------------------------------------------------- NCALC = -1 DO 80 I = 2,NB IF (RATIO .GE. C) GO TO 500 BK(I) = RATIO/EX TWONU = TWONU+TWO RATIO = TWONU 80 CONTINUE NCALC = 1 GO TO 420 ELSE C-------------------------------------------------------------------- C 1.0E-10 .LT. X .LE. 1.0 C-------------------------------------------------------------------- C = ONE X2BY4 = EX*EX/FOUR P0 = HALF*P0 Q0 = HALF*Q0 D1 = -ONE D2 = ZERO BK1 = ZERO BK2 = ZERO F1 = F0 F2 = P0 100 D1 = D1+TWO D2 = D2+ONE D3 = D1+D3 C = X2BY4*C/D2 F0 = (D2*F0+P0+Q0)/D3 P0 = P0/(D2-ENU) Q0 = Q0/(D2+ENU) T1 = C*F0 T2 = C*(P0-D2*F0) BK1 = BK1+T1 BK2 = BK2+T2 IF ((ABS(T1/(F1+BK1)) .GT. EPS) .OR. 1 (ABS(T2/(F2+BK2)) .GT. EPS)) GO TO 100 BK1 = F1+BK1 BK2 = TWO*(F2+BK2)/EX IF (IZE .EQ. 2) THEN D1 = EXP(EX) BK1 = BK1*D1 BK2 = BK2*D1 END IF WMINF = ESTF(1)*EX+ESTF(2) END IF ELSE IF (EPS*EX .GT. ONE) THEN C-------------------------------------------------------------------- C X .GT. ONE/EPS C-------------------------------------------------------------------- NCALC = NB BK1 = ONE / (D*SQRT(EX)) DO 110 I = 1, NB BK(I) = BK1 110 CONTINUE GO TO 500 ELSE C-------------------------------------------------------------------- C X .GT. 1.0 C-------------------------------------------------------------------- TWOX = EX+EX BLPHA = ZERO RATIO = ZERO IF (EX .LE. FOUR) THEN C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 .LE. X .LE. 4.0 C-------------------------------------------------------------------- D2 = AINT(ESTM(1)/EX+ESTM(2)) M = INT(D2) D1 = D2+D2 D2 = D2-HALF D2 = D2*D2 DO 120 I = 2,M D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) 120 CONTINUE C-------------------------------------------------------------------- C Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward C recurrence and K(ALPHA,X) from the wronskian C-------------------------------------------------------------------- D2 = AINT(ESTM(3)*EX+ESTM(4)) M = INT(D2) C = ABS(ENU) D3 = C+C D1 = D3-ONE F1 = XMIN F0 = (TWO*(C+D2)/EX+HALF*EX/(C+D2+ONE))*XMIN DO 130 I = 3,M D2 = D2-ONE F2 = (D3+D2+D2)*F0 BLPHA = (ONE+D1/D2)*(F2+BLPHA) F2 = F2/EX+F1 F1 = F0 F0 = F2 130 CONTINUE F1 = (D3+TWO)*F0/EX+F1 D1 = ZERO T1 = ONE DO 140 I = 1,7 D1 = C*D1+P(I) T1 = C*T1+Q(I) 140 CONTINUE P0 = EXP(C*(A+C*(P(8)-C*D1/T1)-LOG(EX)))/EX F2 = (C+HALF-RATIO)*F1/EX BK1 = P0+(D3*F0-F2+F0+BLPHA)/(F2+F1+F0)*P0 IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(3)*EX+ESTF(4) ELSE C-------------------------------------------------------------------- C Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by backward C recurrence, for X .GT. 4.0 C-------------------------------------------------------------------- DM = AINT(ESTM(5)/EX+ESTM(6)) M = INT(DM) D2 = DM-HALF D2 = D2*D2 D1 = DM+DM DO 160 I = 2,M DM = DM-ONE D1 = D1-TWO D2 = D2-D1 RATIO = (D3+D2)/(TWOX+D1-RATIO) BLPHA = (RATIO+RATIO*BLPHA)/DM 160 CONTINUE BK1 = ONE/((D+D*BLPHA)*SQRT(EX)) IF (IZE .EQ. 1) BK1 = BK1*EXP(-EX) WMINF = ESTF(5)*(EX-ABS(EX-ESTF(7)))+ESTF(6) END IF C-------------------------------------------------------------------- C Calculation of K(ALPHA+1,X) from K(ALPHA,X) and C K(ALPHA+1,X)/K(ALPHA,X) C-------------------------------------------------------------------- BK2 = BK1+BK1*(ENU+HALF-RATIO)/EX END IF C-------------------------------------------------------------------- C Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, C K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 C-------------------------------------------------------------------- NCALC = NB BK(1) = BK1 IF (IEND .EQ. 0) GO TO 500 J = 2-K IF (J .GT. 0) BK(J) = BK2 IF (IEND .EQ. 1) GO TO 500 M = MIN(INT(WMINF-ENU),IEND) DO 190 I = 2,M T1 = BK1 BK1 = BK2 TWONU = TWONU+TWO IF (EX .LT. ONE) THEN IF (BK1 .GE. (XINF/TWONU)*EX) GO TO 195 GO TO 187 ELSE IF (BK1/EX .GE. XINF/TWONU) GO TO 195 END IF 187 CONTINUE BK2 = TWONU/EX*BK1+T1 ITEMP = I J = J+1 IF (J .GT. 0) BK(J) = BK2 190 CONTINUE 195 M = ITEMP IF (M .EQ. IEND) GO TO 500 RATIO = BK2/BK1 MPLUS1 = M+1 NCALC = -1 DO 410 I = MPLUS1,IEND TWONU = TWONU+TWO RATIO = TWONU/EX+ONE/RATIO J = J+1 IF (J .GT. 1) THEN BK(J) = RATIO ELSE IF (BK2 .GE. XINF/RATIO) GO TO 500 BK2 = RATIO*BK2 END IF 410 CONTINUE NCALC = MAX(MPLUS1-K,1) IF (NCALC .EQ. 1) BK(1) = BK2 IF (NB .EQ. 1) GO TO 500 420 J = NCALC+1 DO 430 I = J,NB IF (BK(NCALC) .GE. XINF/BK(I)) GO TO 500 BK(I) = BK(NCALC)*BK(I) NCALC = I 430 CONTINUE END IF 500 RETURN C---------- Last line of RKBESL ---------- END DOUBLE PRECISION FUNCTION DGAMMA(X) C---------------------------------------------------------------------- C C This routine calculates the GAMMA function for a real argument X. C Computation is based on an algorithm outlined in reference 1. C The program uses rational functions that approximate the GAMMA C function to at least 20 significant decimal digits. Coefficients C for the approximation over the interval (1,2) are unpublished. C Those for the approximation for X .GE. 12 are from reference 2. C The accuracy achieved depends on the arithmetic system, the C compiler, the intrinsic functions, and proper selection of the C machine-dependent constants. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta - radix for the floating-point representation C maxexp - the smallest positive power of beta that overflows C XBIG - the largest argument for which GAMMA(X) is representable C in the machine, i.e., the solution to the equation C GAMMA(XBIG) = beta**maxexp C XINF - the largest machine representable floating-point number; C approximately beta**maxexp C EPS - the smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMININ - the smallest positive floating-point number such that C 1/XMININ is machine representable C C Approximate values for some important machines are: C C beta maxexp XBIG C C CRAY-1 (S.P.) 2 8191 966.961 C Cyber 180/855 C under NOS (S.P.) 2 1070 177.803 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 35.040 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 171.624 C IBM 3033 (D.P.) 16 63 57.574 C VAX D-Format (D.P.) 2 127 34.844 C VAX G-Format (D.P.) 2 1023 171.489 C C XINF EPS XMININ C C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 C IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 C VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns the value XINF for singularities or C when overflow would occur. The computation is believed C to be free of underflow and overflow. C C C Intrinsic functions required are: C C INT, DBLE, EXP, LOG, REAL, SIN C C C References: "An Overview of Software Development for Special C Functions", W. J. Cody, Lecture Notes in Mathematics, C 506, Numerical Analysis Dundee, 1975, G. A. Watson C (ed.), Springer Verlag, Berlin, 1976. C C Computer Approximations, Hart, Et. Al., Wiley and C sons, New York, 1968. C C Latest modification: October 12, 1989 C C Authors: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,N LOGICAL PARITY DOUBLE PRECISION 1 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, 2 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, 1 SQRTPI/0.9189385332046727417803297D0/, 2 PI/3.1415926535897932384626434D0/ C---------------------------------------------------------------------- C Machine dependent parameters C---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, 1 XINF/1.79D308/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C approximation over (1,2). C---------------------------------------------------------------------- DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ C---------------------------------------------------------------------- C Coefficients for minimax approximation over (12, INF). C---------------------------------------------------------------------- DATA C/-1.910444077728D-03,8.4171387781295D-04, 1 -5.952379913043012D-04,7.93650793500350248D-04, 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, 3 5.7083835261D-03/ C---------------------------------------------------------------------- C Statement functions for conversion between integer and float C---------------------------------------------------------------------- CONV(I) = DBLE(I) PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y .LE. ZERO) THEN C---------------------------------------------------------------------- C Argument is negative C---------------------------------------------------------------------- Y = -X Y1 = AINT(Y) RES = Y - Y1 IF (RES .NE. ZERO) THEN IF (Y1 .NE. AINT(Y1*HALF)*TWO) PARITY = .TRUE. FACT = -PI / SIN(PI*RES) Y = Y + ONE ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Argument is positive C---------------------------------------------------------------------- IF (Y .LT. EPS) THEN C---------------------------------------------------------------------- C Argument .LT. EPS C---------------------------------------------------------------------- IF (Y .GE. XMININ) THEN RES = ONE / Y ELSE RES = XINF GO TO 900 END IF ELSE IF (Y .LT. TWELVE) THEN Y1 = Y IF (Y .LT. ONE) THEN C---------------------------------------------------------------------- C 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- Z = Y Y = Y + ONE ELSE C---------------------------------------------------------------------- C 1.0 .LT. argument .LT. 12.0, reduce argument if necessary C---------------------------------------------------------------------- N = INT(Y) - 1 Y = Y - CONV(N) Z = Y - ONE END IF C---------------------------------------------------------------------- C Evaluate approximation for 1.0 .LT. argument .LT. 2.0 C---------------------------------------------------------------------- XNUM = ZERO XDEN = ONE DO 260 I = 1, 8 XNUM = (XNUM + P(I)) * Z XDEN = XDEN * Z + Q(I) 260 CONTINUE RES = XNUM / XDEN + ONE IF (Y1 .LT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- RES = RES / Y1 ELSE IF (Y1 .GT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 2.0 .LT. argument .LT. 12.0 C---------------------------------------------------------------------- DO 290 I = 1, N RES = RES * Y Y = Y + ONE 290 CONTINUE END IF ELSE C---------------------------------------------------------------------- C Evaluate for argument .GE. 12.0, C---------------------------------------------------------------------- IF (Y .LE. XBIG) THEN YSQ = Y * Y SUM = C(7) DO 350 I = 1, 6 SUM = SUM / YSQ + C(I) 350 CONTINUE SUM = SUM/Y - Y + SQRTPI SUM = SUM + (Y-HALF)*LOG(Y) RES = EXP(SUM) ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Final adjustments and return C---------------------------------------------------------------------- IF (PARITY) RES = -RES IF (FACT .NE. ONE) RES = FACT / RES 900 DGAMMA = RES RETURN C ---------- Last line of GAMMA ---------- END C SUBROUTINE INVERT(H,HI,A,X,N,IN,N1,IFAULT) C C INVERT SYMMETRIC POSITIVE DEFINITE MATRIX C C ALGORITHM 9 OF NASH, "COMPACT NUMERICAL METHODS FOR COMPUTERS", C SECTION 8.2, P.84. C C H(N,N): MATRIX TO BE INVERTED (INPUT) C HI(N,N): INVERSE OF H (OUTPUT) C N: DIMENSION C IN: LEAD DIMENSION OF ARRAYS H AND IH, AT LEAST N C N1: MUST BE DEFINED IN CALLING PROGRAM TO BE N*(N+1)/2 C A: ARRAY OF DIMENSION AT LEAST N1, USED AS WORKSPACE C X: ARRAY OF DIMENSION AT LEAST N, USED AS WORKSPACE C IFAULT: CONTAINS 0 ON SATISFACTORY COMPLETION, 1 OTHERWISE C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION H(IN,IN),HI(IN,IN),A(N1),X(N) DO 10 I=1,N DO 10 J=1,I K=I*(I-1)/2+J 10 A(K)=H(I,J) DO 20 K=1,N S=A(1) IF(S.LE.0)GOTO100 M=1 DO 40 I=2,N IQ=M M=M+I T=A(IQ+1) X(I)=-T/S IF(I.GT.N+1-K)X(I)=-X(I) IF(IQ+2.GT.M)GOTO40 DO 50 J=IQ+2,M A(J-I)=A(J)+T*X(J-IQ) 50 CONTINUE 40 CONTINUE IQ=IQ-1 A(M)=1/S DO 60 I=2,N A(IQ+I)=X(I) 60 CONTINUE 20 CONTINUE C SUCCESSFUL CONCLUSION: WRITE RESULT INTO HI DO 70 I=1,N DO 70 J=1,I K=I*(I-1)/2+J HI(I,J)=A(K) 70 HI(J,I)=A(K) IFAULT=0 RETURN C UNSUCCESSFUL: SET IFAULT=1 AND LET HI BE IDENTITY MARTIX 100 DO 110 I=1,N DO 110 J=1,N HI(I,J)=0 110 IF(I.EQ.J)HI(I,J)=1 IFAULT=1 RETURN END C