Commit b449b746 authored by agebhard's avatar agebhard
Browse files

Error in bk_ponts found: covmat was not generated and all 0 !!!!

parent 5a48436b
......@@ -110,20 +110,12 @@ bk.points <- function(point.obj,
# ans<-krige.solve(s$x,s$y,point.obj$x,point.obj$y,
# at,covmat,c0vec,c0,trend,rsearch,nsmin,nsmax,mode)
retval<-list(x=ans$datall[1:npts],
y=ans$datall[(npts+1):(2*npts)],
z=ans$datall[(2*npts+1):(3*npts)],
var=ans$datall[(3*npts+1):(4*npts)])
retval$z[ans$napts==1]<-NA
retval$var[ans$napts==1]<-NA
retval
retval<-list(x=ans$xp,
y=ans$yp,
z=ans$zp,
var=ans$varp,
done=ans$dop,
snb=matrix(ans$bits[(npts+2):(nz+n*npts+1)],nrow=n,ncol=npts,byrow=F),
# snb=matrix(ans$bits[(npts+2):(nz+n*npts+1)],nrow=n,ncol=npts,byrow=F),
lambda=matrix(ans$lambda,nrow=n,ncol=npts,byrow=FALSE),
lambda0=ans$lambd0
)
......@@ -132,5 +124,6 @@ bk.points <- function(point.obj,
retval$lambda0[retval$done<=0] <- NA
retval$data<-point.obj
retval$at<-at
retval
}
......@@ -104,7 +104,7 @@ slatec/src/xerprn.o \
slatec/src/xermsg.o
OBJS=$(LINPACK) $(BLAS) $(LAPACK) $(SLATEC) \
#OBJS=$(LINPACK) $(BLAS) $(LAPACK) $(SLATEC) \
OBJS=$(LINPACK) $(LAPACK) $(SLATEC) \
bkwrp.o \
......
......@@ -48,6 +48,8 @@ void bk_point(double *xp,
cov0;
int *indsnb, *indsnw, *indsrt, *ipiv, *ipvt, *iwork, *inddop, *nap, *snb;
int i,j;
double dst;
usesbbt=searchnb[0];
snb=searchnb++;
......@@ -119,6 +121,15 @@ void bk_point(double *xp,
mu =(double *) R_alloc((*ntrend),sizeof(double));
#endif
/* populate covariance matrix */
if(*covtype!=0){
for(i=0;i<(*n);i++)
for(j=i;j<(*n);j++){
dst = sqrt((lon[i]-lon[j])*(lon[i]-lon[j])+(lat[i]-lat[j])*(lat[i]-lat[j]));
covmat[i+j*(*n)]=F77_CALL(covfn)(covtype,covpar,&dst);
}
} /* else { error !! } */
F77_CALL(bk)(xp,
yp,
dop,
......
......@@ -38,6 +38,7 @@ void bk_point(double *xp,
int *ierr,
int *glsmth);
double F77_NAME(covfn)(int *covtype, double *covpar, double *dst);
void F77_NAME(bk)(double *xp,
double *yp,
......
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