code 34441;
procedure GSSNEWTON(M, N, PAR, RV, JJINV, FUNCT, JACOBIAN,
IN, OUT);
value M, N; integer M, N;
array PAR, RV, JJINV, IN, OUT;
boolean procedure FUNCT;
procedure JACOBIAN;
begin integer I, J, INR, MIT, TEXT,
IT, ITMAX, INRMAX, TIM, FEVAL, FEVALMAX;
real RHO, RES1, RES2, RN, RELTOLPAR, ABSTOLPAR, ABSTOLRES,
STAP, NORMX;
boolean CONV, TESTTHF, DAMPING ON;
array JAC[1:M + 1,1:N], PR, AID, SOL[1 : N], FU2[1 : M],
AUX[2 : 5];
integer array CI[1:N];
boolean procedure LOC FUNCT(M, N, PAR, RV);
value M, N; integer M, N; array PAR, RV;
begin LOC FUNCT:= TEST THF:= FUNCT(M, N, PAR, RV)
and TEST THF; FEVAL:= FEVAL + 1
end LOC FUNCT;
ITMAX:= FEVALMAX:= IN[5]; AUX[2]:= N * IN[0]; TIM:= IN[7];
RELTOLPAR:= IN[1] ** 2; ABSTOLPAR:= IN[2] ** 2;
ABSTOLRES:= IN[4] ** 2; INRMAX:= IN[6];
DUPVEC(1, N, 0, PR, PAR);
if M < N then
for I:= 1 step 1 until N do JAC[M + 1, I]:= 0;
TEXT:= 4; MIT:= 0; TEST THF:= true;
RES2:= STAP:= OUT[5]:= OUT[6]:= OUT[7]:= 0;
FUNCT(M, N, PAR, FU2); RN:= VECVEC(1, M, 0, FU2, FU2);
OUT[3]:= SQRT(RN); FEVAL:= 1; DAMPING ON:= false;
for IT:= 1, IT + 1 while IT <= ITMAX and
FEVAL < FEVALMAX do
begin OUT[5]:= IT; JACOBIAN(M, N, PAR, FU2, JAC, LOCFUNCT);
if not TEST THF then
begin TEXT:= 7; "GO TO" FAIL end;
LSQORTDEC(JAC, M, N, AUX, AID, CI);
if AUX[3] ^= N then
begin TEXT:= 5; "GO TO" FAIL end;
LSQSOL(JAC, M, N, AID, CI, FU2); DUPVEC(1, N, 0, SOL, FU2);
STAP:= VECVEC(1, N, 0, SOL, SOL);
RHO:= 2; NORMX:= VECVEC(1, N, 0, PAR, PAR);
if STAP > RELTOLPAR * NORMX + ABSTOLPAR
or IT = 1 and STAP > 0 then
begin for INR:= 0, INR + 1
while if INR = 1 then DAMPING ON or RES2 >= RN
else not CONV and (RN <= RES1 or RES2 < RES1) do
begin comment DAMPING STOPS WHEN
R0 > R1 "AND" R1 <= R2 (BEST RESULT IS X1, R1)
WITH X1 = X0 + I * DX, I:= 1, .5, .25, .125, ETC. ;
RHO:= RHO / 2; if INR > 0 then
begin RES1:= RES2; DUPVEC(1, M, 0, RV, FU2);
DAMPING ON:= INR > 1
end;
for I:= 1 step 1 until N do
PR[I]:= PAR[I] - SOL[I] * RHO;
FEVAL:= FEVAL + 1;
if not FUNCT(M, N, PR, FU2) then
begin TEXT:= 6; "GO TO" FAIL end;
RES2:= VECVEC(1, M, 0, FU2, FU2); CONV:= INR >= INRMAX
end DAMPING OF STEP VECTOR;
if CONV then
begin comment RESIDUE CONSTANT; MIT:= MIT + 1;
if MIT < TIM then CONV:= false
end else MIT:= 0;
if INR > 1 then
begin RHO:= RHO * 2; ELMVEC(1, N, 0, PAR, SOL, - RHO);
RN:= RES1; if INR > 2 then OUT[7]:= IT
end else
begin DUPVEC(1, N, 0, PAR, PR); RN:= RES2;
DUPVEC(1, M, 0, RV, FU2)
end;
if RN <= ABSTOLRES then
begin TEXT:= 1; ITMAX:= IT end else
if CONV and INRMAX > 0 then
begin TEXT:= 3; ITMAX:= IT end
else DUPVEC(1, M, 0, FU2, RV)
end ITERATION WITH DAMPING AND TESTS else
begin TEXT:= 2; RHO:= 1; ITMAX:= IT end
end OF ITERATIONS;
LSQINV(JAC, N, AID, CI);
for I:= 1 step 1 until N do
begin JJINV[I,I]:= JAC[I,I];
for J:= I + 1 step 1 until N do
JJINV[I,J]:= JJINV[J,I]:= JAC[I,J]
end CALCULATION OF INVERSE MATRIX OF NORMAL EQUATIONS;
FAIL :
OUT[6]:= SQRT(STAP) * RHO; OUT[2]:= SQRT(RN); OUT[4]:= FEVAL;
OUT[1]:= TEXT; OUT[8]:= AUX[3]; OUT[9]:= AUX[5]
end GSSNEWTON;
eop