diff --git a/R/bk.grid.R b/R/bk.grid.R index b564ed44ee439a21850f9618711a39909a1d83db..1c92408eab41ca1ac64dc0f669d69a6c1fbb276d 100644 --- a/R/bk.grid.R +++ b/R/bk.grid.R @@ -300,7 +300,7 @@ bk.grid <- function(point.obj, var.mod.obj, xsw=NULL,ysw=NULL,xne=NULL,yne=NULL, dx=NULL,dy=NULL,nx=NULL,ny=NULL, - angle=NULL, + angle=0, maxdist = NULL, extrap = F, border=NULL, @@ -360,8 +360,8 @@ bk.grid <- function(point.obj, method <-"gqr" } method<-switch(method,direct=2,gqr=1,ols=0) - lwork <- glsfit.workquery(n,ntrend,method) - +# lwork <- glsfit.workquery(n,ntrend,method) + lwork <- 3000 if(prior$ntr!=ntrend) stop("model order of priors does not match!") npr<-prior$n diff --git a/makefile b/makefile index 66c01799146f3fc29b5fa55d260fc454f289512f..c48fe26debd695f0400a432b1bfa335a0d62b34b 100644 --- a/makefile +++ b/makefile @@ -1,6 +1,6 @@ all: - (cd .. ; R INSTALL -l $(R_LIBS) --create-lib-so baykrig) + (cd ~/ ; R INSTALL -l $(R_LIBS) baykrig) clean: rm -f src/*.o src/*/*.o src/*/*/*.o src/*.so diff --git a/man/bk.grid.Rd b/man/bk.grid.Rd index 2c7742aa94c1757758d643f2a2c468a3645f2ae6..566726baf68dca6c62fb6d21d1752e1dd735ccfb 100644 --- a/man/bk.grid.Rd +++ b/man/bk.grid.Rd @@ -57,33 +57,32 @@ bk.grid(point.obj, at, prior, var.mod.obj, xsw=NULL, ysw=NULL, xne=NULL, yne=NUL \examples{ \testonly{ # prepare variables from other example pages if they are not already there: - require(rgeostat) - if(length(ls(pat="leman.78"))==0)data(leman) - if(length(ls(pat="leman.bank"))==0)data(leman.bank) - if(length(ls(pat="leman.pt"))==0) leman.pt<-point(leman.88) - if(length(ls(pat="leman.pr"))==0) leman.pr<-pair(leman.pt,maxdist=30) - if(length(ls(pat="leman.ev"))==0) - leman.ev<-est.variogram(leman.pt,leman.pr,"cadpbm") - if(length(ls(pat="leman.sph"))==0) - leman.sph<-fit.variogram("spherical",leman.ev,0.1,0.2,20,iter=50,plot.it=T) - if(length(ls(pat="leman.prior"))==0) example(empirical.prior) + if(length(ls(pat="leman.prior.1"))==0) example(empirical.prior) } -leman.bk <- - bk.grid(point = leman.pt, at = "cadpbm", - prior=leman.prior,var.mod.obj = leman.sph, +leman.bk1 <- + bk.grid(point = leman.88.pt, at = "cadpbm", + prior=leman.prior.1,var.mod.obj = leman.88.vmsph, xsw=min(leman.bank$x),ysw=min(leman.bank$y), xne=max(leman.bank$x), yne=max(leman.bank$y), nx=100, ny=100, trend=1, rsearch = 10, extrap = F,border=leman.bank, duplicate="mean") -plot(leman.bk) +plot(leman.bk1) + +leman.bk0 <- + bk.grid(point = leman.88.pt, at = "cadpbm", + prior=leman.prior.0,var.mod.obj = leman.88.vmsph, + xsw=min(leman.bank$x),ysw=min(leman.bank$y), + xne=max(leman.bank$x), yne=max(leman.bank$y), + nx=100, ny=100, trend=0, rsearch = 10, + extrap = F,border=leman.bank, duplicate="mean") +plot(leman.bk0) + # compare with require(rgeostat) -leman.krige <- - krige.grid(point = leman.pt, at = "cadpbm", - var.mod.obj = leman.sph, - xsw=min(leman.bank$x),ysw=min(leman.bank$y), - xne=max(leman.bank$x), yne=max(leman.bank$y), - nx=100, ny=100, trend=1, rsearch = 10, - extrap = F, border=leman.bank) +leman.krg<-krige.grid.c(min(leman.bank$x),min(leman.bank$y), + max(leman.bank$x), max(leman.bank$y), nx=100, ny=100, point = + leman.88.pt, at = "cadpbm", var.mod.obj = leman.88.vmsph, + rsearch = 7.5,extrap = F,border=leman.bank) + plot(leman.krg,show.snb=T) } \keyword{ baykrig } diff --git a/man/empirical.prior.Rd b/man/empirical.prior.Rd index aec76c2fb80a83540c001eb5ee64b4f6aa3c3ca9..03934889e7b0e99d5c79507dfa9dda85929435af 100644 --- a/man/empirical.prior.Rd +++ b/man/empirical.prior.Rd @@ -37,23 +37,49 @@ empirical.prior(x, formula=NULL, var.mod=NULL, prior=NULL, method="gqr", namx=NU \examples{ \testonly{ # prepare variables from other example pages if they are not already there: +} +# need some preparation steps using sgeostat functions require(rgeostat) +# the data set if(length(ls(pat="leman.78"))==0)data(leman) - if(length(ls(pat="leman.bank"))==0) data(leman.bank) - if(length(ls(pat="leman.pt"))==0) leman.pt<-point(leman.88) - if(length(ls(pat="leman.pr"))==0) leman.pr<-pair(leman.pt,maxdist=30) - if(length(ls(pat="leman.ev"))==0) - leman.ev<-est.variogram(leman.pt,leman.pr,"cadpbm") - if(length(ls(pat="leman.sph"))==0) - leman.sph<-fit.variogram("spherical",leman.ev,0.1,0.2,20,iter=50,plot.it=T) -} + if(length(ls(pat="leman.bank"))==0)data(leman.bank) +# data(leman) contains three measurement grids: +# 1978, 1983, 1988 +# generate the point objects: + if(length(ls(pat="leman.78.pt"))==0) leman.78.pt<-point(leman.78) + if(length(ls(pat="leman.83.pt"))==0) leman.83.pt<-point(leman.83) + if(length(ls(pat="leman.88.pt"))==0) leman.88.pt<-point(leman.88) +# use 1988 data for variogram estimation, star with pairs + if(length(ls(pat="leman.88.pr.30"))==0) + leman.88.pr.30<-pair(leman.88.pt,maxdist=30) +# do the estimation + if(length(ls(pat="leman.88.ev30"))==0) + leman.88.ev30<-est.variogram(leman.88.pt,leman.88.pr.30,"cadpbm") +# fit a sperical variogram + if(length(ls(pat="leman.88.vmsph"))==0) + leman.88.vmsph<-fit.variogram("spherical",leman.88.ev30,0.1,0.2,25,plot.it=T) +plot(leman.88.ev30,var.mod.obj=leman.88.vmsph) + + + +# now prepare the prior guesses: +# simple linear trend +leman.prior.1<-empirical.prior(leman.83,cadpbm~x+y,leman.88.vmsph,duplicate="mean") + +leman.prior.1<-empirical.prior(leman.78,cadpbm~x+y,leman.88.vmsph,prior=leman.prior.1,duplicate="mean") + +# or constant mean model +leman.prior.0<-empirical.prior(leman.83,cadpbm~1,leman.88.vmsph,duplicate="mean",namx="x",namy="y") + +leman.prior.0<-empirical.prior(leman.78,cadpbm~1,leman.88.vmsph,prior=leman.prior.0,duplicate="mean",namx="x",namy="y") + # (some outliers have to be removed:) -leman.prior<-empirical.prior(x = leman.83[-c(69,209),], - formula = cadpbm ~ x + y + 1, - var.mod=leman.sph, duplicate="mean") -leman.prior<-empirical.prior(x = leman.78[-c(69,208),], - formula = cadpbm ~ x + y + 1, - var.mod=leman.sph, prior=leman.prior, - duplicate="mean") +#leman.prior<-empirical.prior(x = leman.83[-c(69,209),], +# formula = cadpbm ~ x + y + 1, +# var.mod=leman.sph, duplicate="mean") +#leman.prior<-empirical.prior(x = leman.78[-c(69,208),], +# formula = cadpbm ~ x + y + 1, +# var.mod=leman.sph, prior=leman.prior, +# duplicate="mean") } \keyword{ baykrig } diff --git a/src/Makevars b/src/Makevars index 28ac24875ea45a599804b35783fb74c070ca0148..ff1781526f04271ef406fc733bdf5ef2a0ad9f89 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,9 +1,20 @@ -OBJS=\ -linpack/dgefa.o \ -linpack/dgedi.o \ -bkwrp.o \ -lsfit.o \ -lapack/util/ilaenv.o \ + +BLAS=blas/dger.o \ +blas/dtrsv.o \ +blas/dtpmv.o \ +blas/dgemm.o \ +blas/dtrsm.o \ +blas/dspr.o \ +blas/dtrmv.o \ +blas/dsymv.o \ +blas/dtrmm.o \ +blas/dsyrk.o \ +blas/dgemv.o \ +blas/dsyr.o \ +blas/dsyr2.o \ +blas/dsyr2k.o \ + +LAPACK=lapack/util/ilaenv.o \ lapack/util/ieeeck.o \ lapack/util/lsame.o \ lapack/double/dpotrf.o \ @@ -76,12 +87,12 @@ lapack/double/dsterf.o \ lapack/double/dsyev.o \ lapack/double/dsytd2.o \ lapack/double/dsytrd.o \ -design.o \ -bk.o \ -bktile.o \ -dsysvq.o \ -dgelse.o \ -slatec/src/j4save.o \ +preload/xerbla.o + +LINPACK= linpack/dgefa.o \ +linpack/dgedi.o + +SLATEC=slatec/src/j4save.o \ slatec/src/fdump.o \ slatec/src/i1mach.o \ slatec/src/xercnt.o \ @@ -90,30 +101,25 @@ slatec/src/dpsort.o \ slatec/src/xerhlt.o \ slatec/src/xersve.o \ slatec/src/xerprn.o \ -slatec/src/xermsg.o \ +slatec/src/xermsg.o + + +OBJS=$(LINPACK) $(BLAS) $(LAPACK) $(SLATEC) \ +bkwrp.o \ +lsfit.o \ +design.o \ +bk.o \ +dsysvq.o \ +dgelse.o \ bkgrid.o \ -preload/xerbla.o \ errmsg.o \ srchnb.o \ dgggle.o \ tools.o \ bkpts.o \ -blas/dger.o \ -blas/dtrsv.o \ -blas/dtpmv.o \ -blas/dgemm.o \ -blas/dtrsm.o \ -blas/dspr.o \ -blas/dtrmv.o \ -blas/dsymv.o \ -blas/dtrmm.o \ -blas/dsyrk.o \ -blas/dgemv.o \ -blas/dsyr.o \ -blas/dsyr2.o \ -blas/dsyr2k.o \ glsfit.o \ covfn.o \ matpr.o \ bk_grid.o +# bktile.o \ diff --git a/src/bk.f b/src/bk.f index 529f2fcf90051634c4512490d7e03e0bd7f516d1..e7be71549ce811fdb1b81593fa87a8e39c2e480a 100644 --- a/src/bk.f +++ b/src/bk.f @@ -42,7 +42,7 @@ c . LDPRIV, . NSMIN, . NSMAX, . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . LDF0WK, @@ -90,7 +90,7 @@ c . LDPRIV, DOUBLE PRECISION LAT0(*),LON0(*),LAT(*),LON(*),Z(*), . COVMAT(LDCOV,*),C0VEC(LDC0,*),COV0,COVPAR(3), . FWORK(LDFWRK,*),F0WORK(LDF0WK,*),MU(LDMU,*), - . Z0(*),LAMBDA(LDLMBD,*),VAR0(*),FWORK2(LDFWRK,*), + . Z0(*),LAMBDA(LDLMBD,*),VAR0(*),FWRK2(LDFWRK,*), . RSEARCH,KWORK(LDKWRK,*),RHSWORK(LDKWRK,*), . FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*), . FPF0WRK(LDFWRK,*),DIST(*),MUPR(LDMPR,*), @@ -255,8 +255,8 @@ c determine search neighbourhood c prepare the design matrix CALL DESIGN(LON,LAT,N,INDSNB,NS,LON0,LAT0,N0,INDDO,NDO, . FWORK,LDFWRK,F0WORK,LDF0WK,NTREND,TREND,IERR) -c make a copy in fwork2 - CALL DSUBMM(FWORK,NS,NTREND,1,1,NS,NTREND,LDFWRK,FWORK2,LDFWRK,1) +c make a copy in fwrk2 + CALL DSUBMM(FWORK,NS,NTREND,1,1,NS,NTREND,LDFWRK,FWRK2,LDFWRK,1) c extract appropriate parts of covariance matrix and Z DO 10 I=1,NS ZSRNB(I)=Z(INDSNB(I)) @@ -270,7 +270,7 @@ c do (generalised) least squares fit in search neighbourhood c CALL GLSFIT(FWORK,N,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV,BETA,U, c . COVBTA,CHLUP,CMINV,WORK1,WORK,LWORK,IPVT,IWORK,IERR) IF (GLSMTH.EQ.0) THEN - CALL LSFIT(FWORK,FWORK2,NS,NTREND,LDFWRK,ZSRNB,BETA, + CALL LSFIT(FWORK,FWRK2,NS,NTREND,LDFWRK,ZSRNB,BETA, . ERRBTA,DEV,ERRDEV, . COVBTA,LDCVBT,SGSQR, . CMINV,LDCINV, @@ -284,7 +284,7 @@ c . COVBTA,CHLUP,CMINV,WORK1,WORK,LWORK,IPVT,IWORK,IERR) RETURN END IF ELSE - CALL GLSFIT(FWORK,FWORK2,NS,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV, + CALL GLSFIT(FWORK,FWRK2,NS,NTREND,LDFWRK,ZSRNB,CVSRNB,LDCOV, . BETA,ERRBTA,DEV,ERRDEV, . COVBTA,LDCVBT,SGSQR, . CHLUP,LDCLUP,CMINV,LDCINV, diff --git a/src/bk_grid.c b/src/bk_grid.c index d53558579569882511dad5f8cbce632c289f3585..0775ddde736693e2316813e88ee0c87d880e581e 100644 --- a/src/bk_grid.c +++ b/src/bk_grid.c @@ -46,14 +46,14 @@ void bk_grid(double *xsw, ldclup=(*n),ldcinv=(*n),ldzg=(*nx); double *covmat, *c0vec, *muwrk, *phiwrk, *beta, errbta, *dev, errdev, *covbta, *cvsrnb, *zsrnb, - *fwork, *fwork2, *f0work, *dist, *kwork, + *fwork, *fwrk2, *f0work, *dist, *kwork, *rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup, *cminv, *work, *ferr, *berr, *mu, *lambda, cov0, lambd0; int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork; -#ifdef TRANSIENT +#ifndef TRANSIENT covmat =Calloc((size_t)(*n)*(*n),double); c0vec =Calloc((size_t)(*n),double); muwrk =Calloc((size_t)((*ntrend)*(*npr)),double); @@ -64,7 +64,7 @@ void bk_grid(double *xsw, cvsrnb =Calloc((size_t)(*n)*(*n),double); zsrnb =Calloc((size_t)(*n),double); fwork =Calloc((size_t)(*n)*(*ntrend),double); - fwork2 =Calloc((size_t)(*n)*(*ntrend),double); + fwrk2 =Calloc((size_t)(*n)*(*ntrend),double); f0work =Calloc((size_t)(*ntrend),double); dist =Calloc((size_t)(*n),double); indsnb =Calloc((size_t)(*n),int); @@ -96,7 +96,7 @@ void bk_grid(double *xsw, cvsrnb =(double *) R_alloc((*n)*(*n),sizeof(double)); zsrnb =(double *) R_alloc((*n),sizeof(double)); fwork =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); - fwork2 =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); + fwrk2 =(double *) R_alloc((*n)*(*ntrend),sizeof(double)); f0work =(double *) R_alloc((*ntrend),sizeof(double)); dist =(double *) R_alloc((*n),sizeof(double)); indsnb =(int *) R_alloc((*n),sizeof(int)); @@ -172,7 +172,7 @@ void bk_grid(double *xsw, nsmin, nsmax, fwork, - fwork2, + fwrk2, &ldfwrk, f0work, dist, @@ -204,8 +204,7 @@ void bk_grid(double *xsw, ierr, glsmth); - -#ifdef TRANSIENT +#ifndef TRANSIENT Free(lambda); Free(mu); Free(iwork); @@ -226,17 +225,55 @@ void bk_grid(double *xsw, Free(indsnb); Free(dist); Free(f0work); - Free(fwork2); + + // Free(fwrk2); Free(fwork); Free(zsrnb); Free(cvsrnb); Free(dev); Free(beta); - Free(covbta); + Free(covbta); // crash bei trend=1 Free(phiwrk); Free(muwrk); Free(c0vec); Free(covmat); + + Free(fwrk2); // crash bei trend=0 + +#else + /* + free(lambda); + free(mu); + free(iwork); + free(berr); + free(ferr); + free(ipiv); + free(ipvt); + free(work); + free(cminv); + free(chlup); + free(fpf0wrk); + free(fpfwork); + free(fpwork); + free(rhswork); + free(kwork); + free(indsrt); + free(indsnw); + free(indsnb); + free(dist); + free(f0work); + free(fwrk2); + free(fwork); + free(zsrnb); + free(cvsrnb); + free(dev); + free(beta); + free(covbta); + free(phiwrk); + free(muwrk); + free(c0vec); + free(covmat); + */ #endif } diff --git a/src/bk_grid.h b/src/bk_grid.h index 5f0425d0b92b1c3b1b004850979082144ad0a860..a3e46a00e4db60042eb0576a6fa265d26255cd7b 100644 --- a/src/bk_grid.h +++ b/src/bk_grid.h @@ -97,7 +97,7 @@ void F77_NAME(bkgrid)(double *xsw, int *nsmin, int *nsmax, double *fwork, - double *fwork2, + double *fwrk2, int *ldfwrk, double *f0work, double *dist, diff --git a/src/bkgrid.f b/src/bkgrid.f index 15aee5beb55c2a81f579ab306e24b6361a11797c..f2c46472ec6c0c919bd9b3a841501d458b9ab23c 100644 --- a/src/bkgrid.f +++ b/src/bkgrid.f @@ -25,7 +25,7 @@ . CVSRNB, . RSEARCH, . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . KWORK, @@ -52,7 +52,7 @@ . COVMAT(LDCOV,*),C0VEC(*),COV0, . RSEARCH,FWORK(LDFWRK,*),F0WORK(*), . KWORK(LDKWRK,*),RHSWORK(*),MU(*), - . LAMBDA(*),FWORK2(LDFWRK,*), + . LAMBDA(*),FWRK2(LDFWRK,*), . COVPAR(*), . FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*), . FPF0WRK(LDFWRK,*),MUPR(LDMPR,*), @@ -143,7 +143,7 @@ c wrapped parameters orig. params(. . INTVEC(11+INTVEC(9)), NSMIN . INTVEC(12+INTVEC(9)), NSMAX . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . DBLVEC(INTVEC(1)+INTVEC(2)+2*INTVEC(8)+2*INTVEC(5)+3),DIST @@ -233,7 +233,7 @@ c wrapped parameters orig. params(. . NSMIN, . NSMAX, . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . DIST, @@ -275,7 +275,7 @@ c wrapped parameters orig. params(. . LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0, . RSEARCH,FWORK(LDFWRK,*),F0WORK(*), . KWORK(LDKWRK,*),RHSWORK(*),MU(*), - . LAMBDA(*),FWORK2(LDFWRK,*), + . LAMBDA(*),FWRK2(LDFWRK,*), . COVPAR(*), . FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*), . FPF0WRK(LDFWRK,*),DIST(*),MUPR(LDMPR,*), @@ -487,7 +487,7 @@ c the main work is now done by BK: . NSMIN, . NSMAX, . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . NTREND, @@ -523,7 +523,6 @@ c the main work is now done by BK: . VAR0, . IERR, . GLSMTH) - c name="xg\0" c call matpr(name,xg,nx,1,nx,dbglvl) c name="yg\0" diff --git a/src/bktile.f b/src/bktile.f index 70babcefb7955968663d045142b10a2f00a6a803..defccc7e90f4fd020c2f7acabc0d74bd0dff1342 100644 --- a/src/bktile.f +++ b/src/bktile.f @@ -27,7 +27,7 @@ c . LDPRIV, . CVSRNB, . RSEARCH, . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . KWORK, @@ -54,7 +54,7 @@ c . LDPRIV, . COVMAT(LDCOV,*),C0VEC(*),COV0, . RSEARCH,FWORK(LDFWRK,*),F0WORK(*), . KWORK(LDKWRK,*),RHSWORK(*),MU(*), - . LAMBDA(*),FWORK2(LDFWRK,*), + . LAMBDA(*),FWRK2(LDFWRK,*), . COVPAR(*), . FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*), . FPF0WRK(LDFWRK,*),MUPR(LDMPR,*), @@ -146,7 +146,7 @@ c . LDPRIV, . INTVEC(11+INTVEC(9)), NSMIN . INTVEC(12+INTVEC(9)), NSMAX . FWORK, - . FWORK2, + . FWRK2, . LDFWRK, . F0WORK, . DBLVEC(INTVEC(1)+INTVEC(2)+2*INTVEC(8)+2*INTVEC(5)+3),DIST @@ -223,7 +223,7 @@ c . LDPRIV, . NSMIN, . NSMAX, . FWORK, - . FWORK2, + . FWRK2, . F0WORK, . DIST, . INDSNB, @@ -256,7 +256,7 @@ c . LDPRIV, . KWORK(NKWORK,*),RHSWORK(NKWORK,*),MU(NTREND,*), . LAMBDA(N,*),X0(*),Y0(*),Z0(*),VAR0(*), . XGWORK(NX,*),YGWORK(NX,*),COVPAR(*), - . FWORK2(N,*) + . FWRK2(N,*) c subroutine for kriging prediction on tiles of a grid c @@ -434,7 +434,7 @@ c fill work arrays with current tile data: c the main work now is done be KRIGE: CALL BK(X0,Y0,DO0,INDDO,LIPT,LON,LAT,Z,N, . COVTYPE,COVPAR,COVMAT,C0,COV0,TREND, - . NTREND,RSEARCH,NSEARCH,NSMIN,NSMAX,FWORK,FWORK2, + . NTREND,RSEARCH,NSEARCH,NSMIN,NSMAX,FWORK,FWRK2, . F0WORK, . DIST,INDSNB,INDSNA,INDSRT,KWORK,NKWORK,RHSWORK, . IPIV,MODE,MU,Z0,LAMBDA,VAR0,IERR)