code 34432;
procedure PRAXIS(N, X, FUNCT, IN, OUT);
value N; integer N;
array X, IN, OUT;
real procedure FUNCT;
begin
commentTHIS PROCEDURE MINIMIZES FUNCT(N,X),WITH THE
PRINCIPAL AXIS METHOD (SEE BRENT,R.P, 1973, ALGORITHMS
FOR MINIMIZATION WITHOUT DERIVATIVES,CH.7);
procedure SORT;
begin integer I, J, K; real S;
for I:= 1 step 1 until N - 1 do
begin K:= I; S:= D[I];
for J:= I+1 step 1 until N do if D[J]>S then
begin K:= J; S:= D[J] end;
if K>I then
begin D[K]:= D[I]; D[I]:= S;
for J:= 1 step 1 until N do
begin S:=V[J,I]; V[J,I]:= V[J,K]; V[J,K]:= S
end
end
end
end SORT
procedure MIN(J, NITS, D2, X1, F1, FK); value J, NITS, FK;
integer J, NITS; real D2, X1, F1; boolean FK;
begin
real procedure FLIN(L); value L; real L;
begin integer I; array T[1:N];
if J > 0 then
begin for I:= 1 step 1 until N do
T[I]:= X[I] + L * V[I,J]
end else
begin comment SEARCH ALONG PARABOLIC SPACE CURVE;
QA:= L * (L - QD1) / (QD0 * (QD0 + QD1));
QB:= (L + QD0) * (QD1 - L) /(QD0 * QD1);
QC:= L * (L + QD0) / (QD1 * (QD0 + QD1));
for I:= 1 step 1 until N do
T[I]:= QA * Q0[I] +QB * X[I] + QC * Q1[I]
end;
NF:= NF + 1; FLIN:= FUNCT(N, T)
end FLIN;
integer K; boolean DZ;
real X2, XM, F0, F2, FM, D1, T2, S, SF1, SX1;
SF1:= F1; SX1:= X1;
K:= 0; XM:= 0; F0:= FM:= FX; DZ:= D2 < RELTOL;
S:= SQRT(VECVEC(1,N,0,X,X));
T2:= M4 * SQRT(ABS(FX) / (if DZ then DMIN else D2)
+ S * LDT) + M2 * LDT; S:= S * M4 + ABSTOL;
if DZ and T2 > S then T2:= S;
if T2<SMALL"THEN"T2:= SMALL;
if T2>0.01*H then T2:= 0.01*H;
if FK"AND"F1<=FM then
begin XM:=X1; FM:= F1 end;
if ^ FK"OR"ABS(X1)<T2"THEN"
begin X1:=if X1>0 then T2"ELSE"-T2;
F1:= FLIN(X1)
end;
if F1<= FM"THEN"
begin XM:= X1; FM:= F1 end;
L0: if DZ then
begin commentEVALUATE FLIN AT ANOTHER POINT
AND ESTIMATE THE SECOND DERIVATIVE;
X2:= if F0 < F1 then -X1 else X1 * 2;
F2:= FLIN(X2); if F2 <= FM then
begin XM:= X2; FM:= F2 end;
D2:=(X2*(F1-F0)-X1*(F2-F0))/(X1*X2*(X1-X2))
end;
commentESTIMATE FIRST DERIVATIVE AT 0;
D1:=(F1-F0)/X1-X1*D2; DZ:=true;
X2:= if D2<=SMALL"THEN"
(if D1<0then H"ELSE"-H)
else -0.5*D1/D2;
if ABS(X2)>H"THEN"X2:=if X2>0then H"ELSE"-H;
L1: F2:=FLIN(X2);
if K<NITS"AND"F2>F0"THEN"
begin K:=K+1;
if F0<F1"AND"X1*X2>0then goto L0;
X2:= 0.5*X2; goto L1
end;
NL:= NL+1;
if F2>FM"THEN"X2:=XM"ELSE"FM:=F2;
D2:=if ABS(X2*(X2-X1))>SMALL"THEN"
(X2*(F1-F0)-X1*(FM-F0))/(X1*X2*(X1-X2))
else if K>0then 0else D2;
if D2<=SMALL"THEN"D2:=SMALL;
X1:=X2; FX:=FM;
if SF1<FX"THEN"
begin FX:=SF1; X1:=SX1 end;
if J>0then ELMVECCOL(1,N,J,X,V,X1)
end MIN;
procedure QUAD;
begin integer I; real L, S;
S:= FX; FX:= QF1; QF1:= S; QD1:= 0;
for I:= 1 step 1 until N do
begin S:=X[I]; X[I]:= L:= Q1[I]; Q1[I]:= S;
QD1:= QD1 + (S - L) ** 2
end;
L:=QD1:=SQRT(QD1); S:= 0;
if (QD0*QD1>DWARF)and NL>=3*N*N"THEN"
begin MIN(0,2,S,L,QF1,true );
QA:= L*(L-QD1)/(QD0*(QD0+QD1));
QB:=(L+QD0)*(QD1-L)/(QD0*QD1);
QC:= L*(L+QD0)/(QD1*(QD0+QD1))
end else
begin FX:= QF1; QA:= QB:= 0; QC:= 1 end;
QD0:= QD1;for I:= 1step 1until N"DO"
begin S:=Q0[I]; Q0[I]:=X[I];
X[I]:= QA*S + QB*X[I]+QC*Q1[I]
end
end QUAD;
boolean ILLC;
integer I, J, K, K2, NL, MAXF, NF, KL, KT, KTM;
real S, SL, DN, DMIN, FX, F1, LDS, LDT, SF, DF, QF1, QD0,
QD1, QA, QB, QC, M2, M4, SMALL, VSMALL, LARGE, VLARGE, SCBD,
LDFAC,T2, MACHEPS, RELTOL, ABSTOL, H;
array V[1:N,1:N], D, Y, Z, Q0, Q1[1:N];
MACHEPS:= IN[0]; RELTOL:= IN[1]; ABSTOL:= IN[2]; MAXF:= IN[5];
H:= IN[6]; SCBD:= IN[7]; KTM:= IN[8]; ILLC:= IN[9] < 0;
SMALL:= MACHEPS ** 2; VSMALL:= SMALL ** 2;
LARGE:= 1/SMALL; VLARGE:= 1/VSMALL;
M2:= RELTOL; M4:= SQRT(M2); SETRANDOM(0.5);
LDFAC:= if ILLC then 0.1 else 0.01;
KT:=NL:=0; NF:=1; OUT[3]:= QF1:=FX:=FUNCT(N,X);
ABSTOL:=T2:= SMALL+ABS(ABSTOL); DMIN:= SMALL;
if H<ABSTOL*100then H:=ABSTOL*100; LDT:=H;
INIMAT(1,N,1,N,V,0);
for I:=1step 1until N"DO"V[I,I]:= 1;
D[1]:= QD0:= 0; DUPVEC(1,N,0,Q1,X);
INIVEC(1,N,Q0,0);
commentMAIN LOOP;
L0: SF:=D[1]; D[1]:= S:= 0;
MIN(1,2,D[1],S,FX,false );
if S <= 0 then MULCOL(1, N, 1, 1, V, V, -1);
if SF <= 0.9 * D[1] or 0.9 * SF >= D[1] then
INIVEC(2,N,D,0);
for K:= 2step 1until N"DO"
begin DUPVEC(1,N,0,Y,X); SF:=FX;
ILLC:= ILLC or KT>0;
L1: KL:=K; DF:= 0; if ILLC then
begin commentRANDOM STOP TO GET OFF
RESULTION VALLEY;
for I:= 1 step 1until N"DO"
begin S:=Z[I]:=(0.1*LDT+T2*10**KT)
*(RANDOM-0.5);
ELMVECCOL(1,N,I,X,V,S)
end;
FX:= FUNCT(N,X); NF:= NF+1
end;
for K2:= K step 1 until N do
begin SL:=FX; S:= 0;
MIN (K2, 2, D[K2], S, FX, false );
S:=if ILLC then D[K2] * (S + Z[K2]) ** 2
else SL-FX;if DF<S"THEN"
begin DF:=S;KL:= K2"END";
end;
if ^ILLC and DF < ABS(100 * MACHEPS * FX) then
begin ILLC:= true; goto L1 end;
for K2:= 1step 1until K-1do
begin S:= 0; MIN(K2, 2, D[K2], S, FX, false ) end;
F1:= FX; FX:= SF; LDS:= 0;
for I:= 1 step 1 until N do
begin SL:= X[I]; X[I]:= Y[I]; SL:= Y[I]:= SL - Y[I];
LDS:= LDS + SL * SL
end; LDS:= SQRT(LDS);
if LDS > SMALL then
begin for I:= KL - 1 step -1 until K do
begin for J:= 1 step 1 until N do
V[J, I + 1]:= V[J,I]; D[I + 1]:= D[I]
end;
D[K]:= 0; DUPCOLVEC(1, N, K, V, Y);
MULCOL(1, N, K, K, V, V, 1 / LDS);
MIN(K, 4, D[K], LDS, F1, true ); if LDS <= 0 then
begin LDS:= LDS; MULCOL(1, N, K, K, V, V, -1) end
end;
LDT:= LDFAC * LDT; if LDT < LDS then LDT:= LDS;
T2:= M2 * SQRT(VECVEC(1, N, 0, X, X)) + ABSTOL;
KT:= if LDT > 0.5 * T2 then 0 else KT + 1;
if KT > KTM then begin OUT[1]:= 0; goto L2 end
end;
QUAD;
DN:= 0;for I:= 1step 1until N"DO"
begin D[I]:= 1/SQRT(D[I]);
if DN<D[I]then DN:=D[I]
end;
for J:= 1step 1until N"DO"
begin S:= D[J]/DN; MULCOL(1,N,J,J,V,V,S)end;
if SCBD>1then
begin S:=VLARGE; for I:=1 step 1until N"DO"
begin SL:= Z[I]:= SQRT(MATTAM(1, N, I, I, V, V));
if SL<M4"THEN"Z[I]:= M4;
if S>SL then S:= SL
end;
for I:=1step 1until N"DO"
begin SL:=S/Z[I];Z[I]:= 1/SL;
if Z[I]>SCBD"THEN"
begin SL:=1/SCBD; Z[I]:= SCBD"END";
MULROW(1, N, I, I, V, V, SL)
end
end;
for I:= 1 step 1 until N do
ICHROWCOL(I + 1, N, I, I, V);
begin array A[1:N,1:N], EM[0:7];
EM[0]:= EM[2]:= MACHEPS;
EM[4]:= 10 * N; EM[6]:= VSMALL;
DUPMAT(1, N, 1, N, A, V);
if QRISNGVALDEC(A, N, N, D, V, EM) ^= 0 then
begin OUT[1]:= 2; goto L2 end;
end;
if SCBD>1then
begin for I:=1step 1until N"DO"
MULROW(1,N,I,I,V,V,Z[I]);
for I:= 1step 1until N"DO"
begin S:= SQRT(TAMMAT(1,N,I,I,V,V));
D[I]:= S*D[I]; S:= 1/S;
MULCOL(1,N,I,I,V,V,S)
end
end;
for I:= 1 step 1 until N do
begin S:= DN * D[I];
D[I]:= if S > LARGE then VSMALL else
if S < SMALL then VLARGE else S ** (-2)
end;
SORT;
DMIN:= D[N]; if DMIN < SMALL then DMIN:= SMALL;
ILLC:= (M2 * D[1]) > DMIN;
if NF < MAXF then goto L0 else OUT[1]:= 1;
L2: OUT[2]:= FX;
OUT[4]:= NF; OUT[5]:= NL; OUT[6]:= LDT
end PRAXIS;
eop