Commit e24b44b3 authored by agebhard's avatar agebhard

covmat was only correct in upper triangle!

parent b449b746
......@@ -252,6 +252,7 @@ c determine search neighbourhood
CALL SRCHNB(LON0,LAT0,DO0,INDDO,N0,NDO,LON,LAT,N,NS,
. RSEARCH,NSEARCH,NSMIN,NSMAX,
. INDSRT,INDSNB,INDSNW,DIST)
write(*,*)"SEARCHNB is:",NS
if (usesbbt.ne.0) then
do 1212 i=1,n
if (indsnb(i).ne.0) then
......@@ -415,10 +416,10 @@ c prepare the right hand side(s)
c solve the system
c LAMBDA = (C0VEC+FWORK*PHIPR*F0WORK)*(COVMAT+FWORK*PHIPR*FWORK')^-1
c name="kwork\0"
c call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
c name="rhswork\0"
c call matpr(name,rhswork,ns,NDO,LDkWrk,dbglvl)
name="kwork\0"
call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
name="rhswork\0"
call matpr(name,rhswork,ns,NDO,LDkWrk,dbglvl)
C BETTER USE DSYSVX ?
CALL DGESV(NS,NDO,KWORK,LDKWRK,IPIV,RHSWORK,LDKWRK,INFO)
......@@ -447,6 +448,7 @@ c extract optimal weights
DO 50 I=1,NS
c write (*,*)i,indsnb(i),j
LAMBDA(INDSNB(I),INDDO(J))=RHSWORK(I,J)
c write (*,*)LAMBDA(INDSNB(I),INDDO(J))
50 CONTINUE
400 CONTINUE
......
......@@ -127,6 +127,8 @@ void bk_point(double *xp,
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);
covmat[j+i*(*n)]=covmat[i+j*(*n)];
}
} /* else { error !! } */
......
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