code 34444;
procedure PEIDE(N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,DERIV,JAC DFDY,
JAC DFDP, CALL YSTART,DATA,MONITOR);
value N,M,NOBS; integer N,M,NOBS,NBP;
array PAR,RES,JTJINV,IN,OUT;
integer array BP;
procedure CALL YSTART,DATA,MONITOR;
boolean procedure DERIV,JAC DFDY,JACDFDP;
begin integer I,J,EXTRA,WEIGHT,NCOL,NROW,AWAY,NPAR,II,JJ,MAX,
NFE,NIS;
real EPS,EPS1,XEND,C,X,T,HMIN,HMAX,RES1,IN3,IN4,FAC3,FAC4;
array AUX[1:3],OBS[1:NOBS],SAVE[-38:6*N],TOBS[0:NOBS],
YP[1:NBP+NOBS,1:NBP+M],YMAX[1:N],Y[1:6*N*(NBP+M+1)],FY[1:N,1:N],
FP[1:N,1:M+NBP];
integer array COBS[1:NOBS];
boolean FIRST,SEC,CLEAN;
real procedure INTERPOL(STARTINDEX,JUMP,K,TOBSDIF);
value STARTINDEX,JUMP,K,TOBSDIF;
integer STARTINDEX,JUMP,K; real TOBSDIF;
begin integer I; real S,R; S:=Y[STARTINDEX]; R:=TOBSDIF;
for I:=1 step 1 until K do
begin STARTINDEX:=STARTINDEX+JUMP;
S:=S+Y[STARTINDEX]*R; R:=R*TOBSDIF
end; INTERPOL:=S
end INTERPOL;
procedure JAC DYDP(NROW,NCOL,PAR,RES,JAC,LOCFUNCT);
value NROW,NCOL; integer NROW,NCOL;
array PAR,RES,JAC; procedure LOCFUNCT;
begin
DUPMAT(1,NROW,1,NCOL,JAC,YP)
end JACOBIAN
boolean procedure FUNCT(NROW,NCOL,PAR,RES);
value NROW,NCOL; integer NROW,NCOL; array PAR,RES;
begin integer L,K,KNEW,FAILS,SAME,KPOLD,N6,NNPAR,J5N,
COBSII;
real XOLD,HOLD,A0,TOLUP,TOL,TOLDWN,TOLCONV,H,CH,CHNEW,
ERROR,DFI,TOBSDIF;
boolean EVALUATE,EVALUATED,DECOMPOSE,CONV;
array A[0:5],DELTA,LAST DELTA,DF,Y0[1:N],JACOB[1:N,1:N];
integer array P[1:N];
real procedure NORM2(AI); real AI;
begin real S,A; S:= "-100;
for I:= 1 step 1 until N do
begin A:= AI/YMAX[I]; S:= S + A * A end;
NORM2:= S
end NORM2;
procedure RESET;
begin if CH < HMIN/HOLD then CH:= HMIN/HOLD else
if CH > HMAX/HOLD then CH:= HMAX/HOLD;
X:= XOLD; H:= HOLD * CH; C:= 1;
for J:= 0 step N until K*N do
begin for I:= 1 step 1 until N do
Y[J+I]:= SAVE[J+I] * C;
C:= C * CH
end;
DECOMPOSE:=true
end RESET;
procedure ORDER;
begin C:= EPS * EPS; J:= (K-1) * (K + 8)/2 - 38;
for I:= 0 step 1 until K do A[I]:= SAVE[I+J];
J:= J + K + 1;
TOLUP := C * SAVE[J];
TOL := C * SAVE[J + 1];
TOLDWN := C * SAVE[J + 2];
TOLCONV:= EPS/(2 * N * (K + 2));
A0:= A[0]; DECOMPOSE:= true;
end ORDER;
procedure EVALUATE JACOBIAN;
begin EVALUATE:= false;
DECOMPOSE:= EVALUATED:= true;
if not JAC DFDY(PAR,Y,X,FY) then
begin SAVE[-3]:=4; goto RETURN end;
end EVALUATE JACOBIAN
procedure DECOMPOSE JACOBIAN;
begin DECOMPOSE:= false;
C:= -A0 * H;
for J:= 1 step 1 until N do
begin for I:= 1 step 1 until N do
JACOB[I,J]:= FY[I,J] * C;
JACOB[J,J]:= JACOB[J,J] + 1
end;
DEC(JACOB,N,AUX,P)
end DECOMPOSE JACOBIAN;
procedure CALCULATE STEP AND ORDER;
begin real A1,A2,A3;
A1:= if K <= 1 then 0 else
0.75 * (TOLDWN/NORM2(Y[K*N+I])) ** (0.5/K);
A2:= 0.80 * (TOL/ERROR) ** (0.5/(K + 1));
A3:= if K >= 5 or FAILS ^= 0
then 0 else
0.70 * (TOLUP/NORM2(DELTA[I] - LAST DELTA[I]))**
(0.5/(K+2));
if A1 > A2 and A1 > A3 then
begin KNEW:= K-1; CHNEW:= A1 end else
if A2 > A3 then
begin KNEW:= K ; CHNEW:= A2 end else
begin KNEW:= K+1; CHNEW:= A3 end
end CALCULATE STEP AND ORDER;
if SEC then begin SEC:=false; goto RETURN end;
NPAR:=M; EXTRA:=NIS:=0; II:=1;
JJ:=if NBP=0 then 0 else 1;
N6:=N*6;
INIVEC(-3,-1,SAVE,0);
INIVEC(N6+1,(6+M)*N,Y,0);
INIMAT(1,NOBS+NBP,1,M+NBP,YP,0);
T:=TOBS[1]; X:=TOBS[0];
CALL YSTART(PAR,Y,YMAX);
HMAX:=TOBS[1]-TOBS[0]; HMIN:=HMAX*IN[1];
EVALUATE JACOBIAN; NNPAR:=N*NPAR;
NEW START:
K:= 1; KPOLD:=0; SAME:= 2; ORDER;
if not DERIV(PAR,Y,X,DF) then
begin SAVE[-3]:=3; goto RETURN end;
H:=SQRT(2 * EPS/SQRT(NORM2 (MATVEC(1,N,I,FY,DF))));
if H > HMAX then H:= HMAX else
if H < HMIN then H:= HMIN;
XOLD:= X; HOLD:= H; CH:= 1;
for I:= 1 step 1 until N do
begin SAVE[I]:=Y[I]; SAVE[N+I]:=Y[N+I]:=DF[I]*H end;
FAILS:= 0;
for L:= 0 while X < XEND do
begin if X + H <= XEND then X:= X + H else
begin H:= XEND-X; X:= XEND; CH:= H/HOLD; C:= 1;
for J:= N step N until K*N do
begin C:= C* CH;
for I:= J+1 step 1 until J+N do
Y[I]:= Y[I] * C
end;
SAME:= if SAME<3 then 3 else SAME+1;
end;
comment PREDICTION;
for L:= 1 step 1 until N do
begin for I:= L step N until (K-1)*N+L do
for J:= (K-1)*N+L step -N until I do
Y[J]:= Y[J] + Y[J+N];
DELTA[L]:= 0
end; EVALUATED:= false;
comment CORRECTION AND ESTIMATION LOCAL ERROR;
for L:= 1,2,3 do
begin if not DERIV(PAR,Y,X,DF) then
begin SAVE[-3]:=3; goto RETURN end;
for I:= 1 step 1 until N do
DF[I]:= DF[I] * H - Y[N+I];
if EVALUATE then EVALUATE JACOBIAN;
if DECOMPOSE then DECOMPOSE JACOBIAN;
SOL(JACOB,N,P,DF);
CONV:= true;
for I:= 1 step 1 until N do
begin DFI:= DF[I];
Y[ I]:= Y[ I] + A0 * DFI;
Y[N+I]:= Y[N+I] + DFI;
DELTA[I]:= DELTA[I] + DFI;
CONV:= CONV and ABS(DFI) < TOLCONV * YMAX[I]
end;
if CONV then
begin ERROR:= NORM2(DELTA[I]);
goto CONVERGENCE
end
end;
comment ACCEPTANCE OR REJECTION;
if not CONV then
begin if not EVALUATED then EVALUATE:= true
else
begin CH:=CH/4; if H<4*HMIN then
begin SAVE[-1]:= SAVE[-1]+10;
HMIN:=HMIN/10;
if SAVE[-1]>40 then goto RETURN
end
end;
RESET
end else CONVERGENCE:
if ERROR > TOL then
begin FAILS:= FAILS + 1;
if H > 1.1 * HMIN then
begin if FAILS > 2 then
begin RESET; goto NEW START
end else
begin CALCULATE STEP AND ORDER;
if KNEW ^= K then
begin K:= KNEW; ORDER end;
CH:= CH * CHNEW; RESET
end
end else
begin if K = 1 then
begin comment VIOLATE EPS CRITERION;
SAVE[-2]:= SAVE[-2] + 1;
SAME:= 4; goto ERROR TEST OK
end;
K:=1; RESET; ORDER; SAME:= 2
end
end else ERROR TEST OK:
begin FAILS:= 0;
for I:= 1 step 1 until N do
begin C:= DELTA[I];
for L:= 2 step 1 until K do
Y[L*N+I]:= Y[L*N+I] + A[L] * C;
if ABS(Y[I]) > YMAX[I] then
YMAX[I]:= ABS(Y[I])
end;
SAME:= SAME-1;
if SAME= 1 then
DUPVEC(1,N,0,LAST DELTA,DELTA) else
if SAME= 0 then
begin CALCULATE STEP AND ORDER;
if CHNEW > 1.1 then
begin
if K ^= KNEW then
begin if KNEW > K then
MULVEC(KNEW*N+1,KNEW*N+N,-KNEW*N,Y,DELTA,
A[K]/KNEW);
K:= KNEW; ORDER
end;
SAME:= K+1;
if CHNEW * H > HMAX
then CHNEW:= HMAX/H;
H:= H * CHNEW; C:= 1;
for J:= N step N until K*N do
begin C:= C * CHNEW;
MULVEC(J+1,J+N,0,Y,Y,C)
end; DECOMPOSE:=true
end
else SAME:= 10
end OF A SINGLE INTEGRATION STEP OF Y;
NIS:=NIS+1;
comment START OF A INTEGRATION STEP OF YP;
if CLEAN then
begin HOLD:=H; XOLD:=X; KPOLD:=K; CH:=1;
DUPVEC(1,K*N+N,0,SAVE,Y)
end else
begin if H^=HOLD then
begin CH:=H/HOLD; C:=1;
for J:=N6+NNPAR step NNPAR until
KPOLD*NNPAR+N6 do
begin C:=C*CH;
for I:=J+1 step 1 until J+NNPAR do
Y[I]:=Y[I]*C
end; HOLD:=H
end;
if K>KPOLD then
INIVEC(N6+K*NNPAR+1,N6+K*NNPAR+NNPAR,Y,0);
XOLD:= X; KPOLD:= K; CH:= 1;
DUPVEC(1,K*N+N,0,SAVE,Y);
EVALUATE JACOBIAN;
DECOMPOSE JACOBIAN;
if not JAC DFDP(PAR,Y,X,FP) then
begin SAVE[-3]:=5; goto RETURN end;
if NPAR>M then INIMAT(1,N,M+1,NPAR,FP,0);
comment PREDICTION;
for L:=0 step 1 until K-1 do
for J:=K-1 step -1 until L do
ELMVEC(J*NNPAR+N6+1,J*NNPAR+N6+NNPAR,NNPAR,Y,Y,1);
comment CORRECTION;
for J:=1 step 1 until NPAR do
begin J5N:=(J+5)*N;
DUPVEC(1,N,J5N,Y0,Y);
for I:=1 step 1 until N do DF[I]:=
H*(FP[I,J]+MATVEC(1,N,I,FY,Y0))
-Y[NNPAR+J5N+I];
SOL(JACOB,N,P,DF);
for L:=0 step 1 until K do
begin I:=L*NNPAR+J5N;
ELMVEC(I+1,I+N,-I,Y,DF,A[L])
end
end
end;
for L:=0 while X>=T do
begin
comment CALCULATION OF A ROW OF THE JACOBIAN
MATRIX AND AN ELEMENT OF THE RESIDUAL
VECTOR;
TOBSDIF:=(TOBS[II]-X)/H; COBSII:=COBS[II];
RES[II]:=INTERPOL(COBSII,N,K,TOBSDIF)-OBS[II];
if not CLEAN then
begin for I:=1 step 1 until NPAR do
YP[II,I]:=INTERPOL(COBSII+(I+5)*N,NNPAR,K,
TOBSDIF);
comment INTRODUCING OF BREAK-POINTS;
if BP[JJ]^=II then else
if FIRST and ABS(RES[II])<EPS1 then
begin NBP:=NBP-1; comment DUPVEC(BP) ;
for I:= JJ step 1 until NBP do
BP[I]:= BP[I + 1]; BP[NBP+1]:=0
end else
begin EXTRA:=EXTRA+1;
if FIRST then PAR[M+JJ]:=OBS[II];
comment INTRODUCING A JACOBIAN ROW AND A
RESIDUAL VECTOR ELEMENT FOR
CONTINUITY REQUIREMENTS;
YP[NOBS+JJ,M+JJ]:=-WEIGHT;
MULROW(1,NPAR,NOBS+JJ,II,YP,YP,WEIGHT);
RES[NOBS+JJ]:=WEIGHT*(RES[II]+OBS[II]-
PAR[M+JJ])
end
end;
if II=NOBS then goto RETURN else
begin T:=TOBS[II+1];
if BP[JJ]=II and JJ<NBP then JJ:=JJ+1;
HMAX:=T-TOBS[II]; HMIN:=HMAX*IN[1]; II:=II+1
end;
end;
comment BREAK-POINTS INTRODUCE NEW INITIAL VALUES
FOR Y AND YP;
if EXTRA>0 then
begin for I:=1 step 1 until N do
begin Y[I]:=INTERPOL(I,N,K,TOBSDIF);
for J:=1 step 1 until NPAR do
Y[I+(J+5)*N]:=INTERPOL(I+(J+5)*N,NNPAR,K,
TOBSDIF)
end;
for L:=1 step 1 until EXTRA do
begin COBSII:=COBS[BP[NPAR-M+L]];
Y[COBSII]:=PAR[NPAR+L];
for I:=1 step 1 until NPAR+EXTRA do
Y[COBSII+(5+I)*N]:=0;
INIVEC(1+NNPAR+(L+5)*N,NNPAR+(L+6)*N,Y,0);
Y[COBSII+(5+NPAR+L)*N]:=1
end;
NPAR:=NPAR+EXTRA; EXTRA:=0;
X:=TOBS[II-1]; EVALUATE JACOBIAN; NNPAR:=N*NPAR;
goto NEW START
end
end
end STEP;
RETURN:
if SAVE[-2]>MAX then MAX:=SAVE[-2];
FUNCT:=SAVE[-1]<=40 and SAVE[-3]=0;
if not FIRST then
MONITOR(1,NCOL,NROW,PAR,RES,WEIGHT,NIS)
end FUNCT;
I:= -39;
for C:= 1,1,9,4,0,2/3,1,1/3,36,20.25,1,6/11,
1,6/11,1/11,84.028,53.778,0.25,.48,1,.7,.2,.02,
156.25, 108.51, .027778, 120/274, 1, 225/274,
85/274, 15/274, 1/274, 0, 187.69, .0047361
do begin I:= I + 1; SAVE[I]:= C end;
DATA(NOBS,TOBS,OBS,COBS); WEIGHT:=1;
FIRST:=SEC:=false; CLEAN:=NBP>0;
AUX[2]:="-12; EPS:=IN[2]; EPS1:="10;
XEND:=TOBS[NOBS]; OUT[1]:=0; BP[0]:=MAX:=0;
comment SMOOTH INTEGRATION WITHOUT BREAK-POINTS;
if not FUNCT(NOBS,M,PAR,RES) then goto ESCAPE;
RES1:=SQRT(VECVEC(1,NOBS,0,RES,RES)); NFE:=1;
if IN[5]=1 then
begin OUT[1]:=1; goto ESCAPE end;
if CLEAN then
begin FIRST:=true; CLEAN:=false;
FAC3:=SQRT(SQRT(IN[3]/RES1)); FAC4:=SQRT(SQRT(IN[4]/RES1));
EPS1:=RES1*FAC4;
if not FUNCT(NOBS,M,PAR,RES) then goto ESCAPE;
FIRST:=false
end else NFE:=0;
NCOL:=M+NBP; NROW:=NOBS+NBP;
SEC:=true;
IN3:=IN[3]; IN4:=IN[4]; IN[3]:=RES1;
begin real W; array AID[1:NCOL,1:NCOL];
WEIGHT:=AWAY:=0;
OUT[4]:=OUT[5]:=W:=0;
for WEIGHT:=(SQRT(WEIGHT)+1)**2 while
WEIGHT^=16 and NBP>0 do
begin if AWAY=0 and W^=0 then
begin comment IF NO BREAK-POINTS WERE OMITTED THEN ONE
FUNCTION EVALUATION IS SAVED;
W:=WEIGHT/W;
for I:=NOBS+1 step 1 until NROW do
begin for J:=1 step 1 until NCOL do
YP[I,J]:=W*YP[I,J];
RES[I]:=W*RES[I]
end; SEC:=true; NFE:=NFE-1
end;
IN[3]:=IN[3]*FAC3*WEIGHT; IN[4]:=EPS1;
MONITOR(2,NCOL,NROW,PAR,RES,WEIGHT,NIS);
MARQUARDT(NROW,NCOL,PAR,RES,AID,FUNCT,JAC DYDP,IN,OUT);
if OUT[1]>0 then goto ESCAPE;
comment THE RELATIVE STARTING VALUE OF LAMBDA IS
ADJUSTED TO THE LAST VALUE OF LAMBDA USED;
AWAY:=OUT[4]-OUT[5]-1;
IN[6]:=IN[6] * 5**AWAY * 2**(AWAY-OUT[5]);
NFE:=NFE+OUT[4];
W:=WEIGHT; EPS1:=(SQRT(WEIGHT)+1)**2*IN[4]*FAC4;
AWAY:=0;
comment USELESS BREAK-POINTS ARE OMITTED;
J:= 0; for J:= J + 1 while J le NBP do
begin if ABS(OBS[BP[J]]+RES[BP[J]]-PAR[J+M])<EPS1
then
begin NBP:=NBP-1; comment DUPVEC (BP) ;
for I:= J step 1 until NBP do
BP[I]:= BP[I + 1];
DUPVEC(J+M,NBP+M,1,PAR,PAR);
J:=J-1; AWAY:=AWAY+1; BP[NBP+1]:=0
end
end;
NCOL:=NCOL-AWAY; NROW:=NROW-AWAY
end;
IN[3]:=IN3; IN[4]:=IN4; NBP:=0; WEIGHT:=1;
MONITOR(2,M,NOBS,PAR,RES,WEIGHT,NIS);
MARQUARDT(NOBS,M,PAR,RES,JTJINV,FUNCT,JAC DYDP,IN,OUT);
NFE:=OUT[4]+NFE
end;
ESCAPE: if OUT[1]=3 then OUT[1]:=2 else
if OUT[1]=4 then OUT[1]:=6;
if SAVE[-3]^=0 then OUT[1]:=SAVE[-3];
OUT[3]:=RES1;
OUT[4]:=NFE;
OUT[5]:=MAX
end PEIDE;
eop