Commit 8fb00b82 authored by agebhard's avatar agebhard

error found! different sorting in vectors leads to DDOT(,)=0 !!!

parent e24b44b3
......@@ -252,7 +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
c write(*,*)"SEARCHNB is:",NS
if (usesbbt.ne.0) then
do 1212 i=1,n
if (indsnb(i).ne.0) then
......@@ -416,9 +416,9 @@ c prepare the right hand side(s)
c solve the system
c LAMBDA = (C0VEC+FWORK*PHIPR*F0WORK)*(COVMAT+FWORK*PHIPR*FWORK')^-1
name="kwork\0"
call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
name="rhswork\0"
c name="kwork\0"
c call matpr(name,kwork,ns,ns,LDkwrk,dbglvl)
c name="rhswork\0"
call matpr(name,rhswork,ns,NDO,LDkWrk,dbglvl)
C BETTER USE DSYSVX ?
......@@ -489,15 +489,17 @@ c C0VEC)
IF (MODE.EQ.2 .OR. MODE.EQ.3) THEN
c work array: RHSWORK = C0VEC + FWORK PHIPR F0WORK
c the same procedure as to prepare the right hand side(s) above:
c (use the same index sorting as in lambda!!)
DO 1300 J=1,NDO
DO 140 I=1,NS
RHSWORK(I,J)=C0VEC(INDSNB(I),inddo(J))+FPF0WRK(I,j)
RHSWORK(indsnb(I),J)=C0VEC(INDSNB(I),inddo(J))+FPF0WRK(I,J)
140 CONTINUE
1300 CONTINUE
DO 500 I=1,NDO
c work array: MU = PHIPR * F0WORK
CALL DGEMV('N',NTREND,NTREND,ONE,PHIWRK,LDPHWK,
. F0WORK(1,I),1,ZERO,MU,1)
. F0WORK(1,I),1,ZERO,MU,1)
c write(*,*)"cov0:",cov0
c name="mu\0"
c call matpr(name,mu,ntrend,1,ldmu,dbglvl)
c name="f0work\0"
......@@ -508,6 +510,8 @@ c name="lambda\0"
c call matpr(name,lambda(1,inddo(i)),n,1,ldlmbd,dbglvl)
c name="c0vec\0"
c call matpr(name,c0vec(1,inddo(i)),n,1,ldc0,dbglvl)
c name="rhswork\0"
c call matpr(name,rhswork(1,inddo(i)),n,1,ldkwrk,dbglvl)
VAR0(INDDO(I)) = COV0 - DDOT(N,LAMBDA(1,INDDO(I)),
. 1,RHSWORK(1,INDDO(I)),1)
......
......@@ -137,9 +137,9 @@ void bk_point(double *xp,
dop,
inddop,
np,
lon,
lat,
z,
lon,
lat,
z,
n,
covtype,
covpar,
......
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