Commit 3bfd84bf authored by alge's avatar alge

still problems with segfqults, but only when closing R

added snb structure. lambda and lambda0 new return values.
new bias=T plot in rgeostat.
plot.krige.map doesn't work for non square grids? for bias plots?
parent 6875d810
......@@ -76,6 +76,8 @@ c . LDPRIV,
. LDLMBD,
. LAMBD0,
. VAR0,
. SNBBIT,
. USESBBT,
. IERR,
. GLSMTH)
......@@ -86,7 +88,7 @@ c . LDPRIV,
. NPR,TYPPR(*),LWORK,IPVT(*),IWORK(*),GLSMTH,
. LDCVBT,LDCLUP,LDCINV,
. LDMPR,LDPHPR,LDPHWK, LDPRIV
. LDF0WK,LDRSWK,LDMU
. LDF0WK,LDRSWK,LDMU,USESBBT,SNBBIT(*)
DOUBLE PRECISION LAT0(*),LON0(*),LAT(*),LON(*),Z(*),
. COVMAT(LDCOV,*),C0VEC(LDC0,*),COV0,COVPAR(3),
. FWORK(LDFWRK,*),F0WORK(LDF0WK,*),MU(LDMU,*),
......@@ -252,6 +254,15 @@ c determine search neighbourhood
CALL SRCHNB(LON0,LAT0,DO0,INDDO,N0,NDO,LON,LAT,N,NS,
. RSEARCH,NSEARCH,NSMIN,NSMAX,
. INDSRT,INDSNB,INDSNW,DIST)
if (usesbbt.ne.0) then
do 1212 i=1,n
if (indsnb(i).ne.0) then
snbbit(indsnb(i))=1
else
snbbit(indsnb(i))=0
end if
1212 continue
end if
c prepare the design matrix
CALL DESIGN(LON,LAT,N,INDSNB,NS,LON0,LAT0,N0,INDDO,NDO,
......
......@@ -37,20 +37,23 @@ void bk_grid(double *xsw,
int *nsmax,
int *lwork,
int *mode,
double *lambda,
double *lambd0,
int *searchnb,
int *ierr,
int *glsmth){
int nz=(*nx)*(*ny), ldcov=(*n), ldc0=(*n),ldphwk=(*ntrend),
ldfwrk=(*n), ldlmbd=(*n), ldkwrk=(*n), // ldcvbt=(*ntrend)*(*ntrend),
ldcvbt=(*ntrend),
ldcvbt=(*ntrend),
ldclup=(*n),ldcinv=(*n),ldzg=(*nx);
double *covmat, *c0vec, *muwrk, *phiwrk, *beta, errbta,
*dev, errdev, *covbta, *cvsrnb, *zsrnb,
*fwork, *fwrk2, *f0work, *dist, *kwork,
*rhswork, *fpwork, *fpfwork, *fpf0wrk, *chlup,
*cminv, *work, *ferr, *berr,
*mu, *lambda, cov0, lambd0;
*mu, // *lambda,
cov0; // lambd0
int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork;
......@@ -86,7 +89,7 @@ void bk_grid(double *xsw,
berr =Calloc((size_t)(*n),double);
iwork =Calloc((size_t)(3*(*n)),int);
mu =Calloc((size_t)(*ntrend),double);
lambda =Calloc((size_t)(*n),double);
// lambda =Calloc((size_t)(*n),double);
#else
covmat =(double *) R_alloc((*n)*(*n),sizeof(double));
c0vec =(double *) R_alloc((*n),sizeof(double));
......@@ -118,7 +121,7 @@ void bk_grid(double *xsw,
berr =(double *) R_alloc((*n),sizeof(double));
iwork =(int *) R_alloc(3*(*n),sizeof(int));
mu =(double *) R_alloc((*ntrend),sizeof(double));
lambda =(double *) R_alloc((*n),sizeof(double));
// lambda =(double *) R_alloc((*n),sizeof(double));
#endif
F77_CALL(bkgrid)(xsw,
......@@ -201,13 +204,14 @@ void bk_grid(double *xsw,
mode,
mu,
lambda,
&lambd0,
&ldlmbd,
lambd0,
searchnb,
ierr,
glsmth);
#ifndef TRANSIENT
Free(lambda);
// Free(lambda);
Free(mu);
Free(iwork);
Free(berr);
......@@ -243,7 +247,7 @@ void bk_grid(double *xsw,
#else
/*
free(lambda);
// free(lambda);
free(mu);
free(iwork);
free(berr);
......
......@@ -40,6 +40,8 @@ void bk_grid(double *xsw,
int *nsmax,
int *lwork,
int *mode,
double *lambda,
double *lambd0,
int *searchnb,
int *ierr,
int *glsmth);
......@@ -124,6 +126,7 @@ void F77_NAME(bkgrid)(double *xsw,
int *mode,
double *mu,
double *lambda,
int *ldlmbd,
double *lambd0,
int *searchnb,
int *ierr,
......
......@@ -78,6 +78,7 @@
. MODE,
. MU,
. LAMBDA,
. LDLMBD,
. LAMBD0,
. BITS,
. IERR,
......@@ -86,18 +87,18 @@
INTEGER NX,NY,NZ,N,COVTYPE,TREND,NTREND,LDZG,
. NSEARCH,NSMIN,NSMAX,MODE,IERR,INDSNB(*),INDSNW(*),
. INDSRT(*),IPIV(*),IPT,EXTRAP,DOG(NX,*),LDKWRK,
. LDCOV,LDFWRK,LDMPR,LDPHPR,LDPHWK,LDCVBT,LWORK,
. LDCOV,LDFWRK,LDMPR,LDPHPR,LDPHWK,LDCVBT,LWORK,LDLMBD,
. LDCLUP,LDCINV,NPR,TYPPR(*),IPVT(*),IWORK(*),GLSMTH,BITS(*)
DOUBLE PRECISION XSW,YSW,XNE,YNE,ANGLE,DX,DY,
. XG(*),YG(*),ZG(LDZG,*),VARG(LDZG,*),
. LON(*),LAT(*),Z(*),COVMAT(LDCOV,*),C0VEC(*),COV0,
. RSEARCH,FWORK(LDFWRK,*),F0WORK(*),
. KWORK(LDKWRK,*),RHSWORK(*),MU(*),
. LAMBDA(*),FWRK2(LDFWRK,*),
. LAMBDA(LDLMBD,*),FWRK2(LDFWRK,*),
. COVPAR(*),
. FPWORK(LDFWRK,*),FPFWORK(LDFWRK,*),
. FPF0WRK(LDFWRK,*),DIST(*),MUPR(LDMPR,*),
. PHIPR(LDPHPR,*),LAMBD0,PHIWRK(LDPHWK,*),
. PHIPR(LDPHPR,*),LAMBD0(LDZG,*),PHIWRK(LDPHWK,*),
. LONPR(*),LATPR(*),MUWRK(*),
. BETA(*),COVBTA(LDCVBT,*),
. WORK(LWORK),
......@@ -152,7 +153,7 @@ c ... ... other work arrays to pass through to KRIGE
c local variables
INTEGER I,J,K,L, INDDO(1), DO0(1), NA0, usesbbt, pcnt
DOUBLE PRECISION DELTA, X0, Y0, Z0, VAR0,DST
DOUBLE PRECISION DELTA, X0, Y0, Z0, VAR0,DST, L0
c debug options
CHARACTER*16 NAME
......@@ -335,10 +336,14 @@ c the main work is now done by BK:
. NTREND,
. Z0,
. NA0,
. LAMBDA,
c . LAMBDA,
. LAMBDA(1,pcnt),
. N,
. LAMBD0,
. L0,
c . LAMBD0(pcnt),
. VAR0,
. BITS(nz+2+(pcnt-1)*n),
. USESBBT,
. IERR,
. GLSMTH)
c name="xg\0"
......@@ -351,10 +356,12 @@ c extract results for this tile
DOG(I,J)=-1
ZG(I,J)=0
VARG(I,J)=0
LAMBD0(I,J)=0
c write(*,*)"x"
ELSE
ZG(I,J)=Z0
VARG(I,J)=VAR0
LAMBD0(I,J)=L0
c write(*,*)"o"
c VARG(I,J)=IWORK(1)*1.0D0
c VARG(I,J)=lambd0
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment