--- /dev/null
+ LONG REAL PROCEDURE ALGOLPRAXIS (LONG REAL VALUE T, MACHEPS, H ;
+ INTEGER VALUE N, PRIN;
+ LONG REAL ARRAY X(*);
+ LONG REAL PROCEDURE F(long real array x(*);
+ integer value n));
+ BEGIN COMMENT:
+ THIS PROCEDURE MINIMIZES THE FONCTION F(X, N) OF N
+ VARIABLES X(1), ... X(N), USING THE PRINCIPAL AXIS METHOD.
+ ON ENTRY X HOLDS A GUESS, ON RETURN IT HOLDS THE ESTIMATED
+ POINT OF MINIMUM, WITH (HOPEFULLY) |ERROR| <
+ SQRT(MACHEPS)*|X| + T, WHERE MACHEPS IS THE MACHINE
+ PRECISION, THE SMALLEST NUMBER SUCH THAT 1 + MACHEPS > 1,
+ T IS A TOLERANCE, AND |.| IS THE 2-NORM. H IS THE MAXIMUM
+ STEP SIZE: SET TO ABOUT THE MAXIMUM EXPECTED DISTANCE FROM
+ THE GUESS TO THE MINIMUM (IF H IS SET TOO SMALL OR TOO
+ LARGE THEN THE INITIAL RATE OF CONVERGENCE WILL BE SLOW).
+ THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS
+ AFTER PROCEDURE QUAD.
+ PRIN CONTROLS THE PRINTING OF INTERMEDIATE RESULTS.
+ IF PRIN = 0, NO RESULTS ARE PRINTED.
+ IF PRIN = 1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR
+ MINIMIZATIONS, AND FINAL X IS PRINTED, BUT INTERMEDIATE
+ X ONLY IF N <= 4.
+ IF PRIN = 2, EIGENVALUES OF A AND SCALE FACTORS ARE ALSO PRINTED.
+ IF PRIN = 3, F AND X ARE PRINTED AFTER EVERY FEW LINEAR MINIMIZATIONS.
+ IF PRIN = 4, EIGENVECTORS ARE ALSO PRINTED.
+ FMIN IS A GLOBAL VARIABLE: SEE PROCEDURE PRINT.
+ RANDOM IS A PARAMETERLESS LONG REAL PROCEDURE WHICH RETURNS
+ A RANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0, 1). ANY
+ INITIALIZATION MUST BE DONE BEFORE THE CALL TO PRAXIS.
+ THE PROCEDURE IS MACHINE-INDEPENDENT, APART FROM THE OUTPUT
+ STATEMENTS AND THE SPECIFICATION OF MACHEPS. WE ASSUME THAT
+ MACHEPS**(—4) DOES NOT OVERFLOW (IF IT DOES THEN MACHEPS MUST
+ BE INCREASED), AND THAT ON FLOATING-POINT UNDERFLOW THE
+ RESULT IS SET TO ZERO;
+
+ LONG REAL PROCEDURE RANDOM(INTEGER VALUE NAUGHT);
+ ALGOL "random";
+
+ PROCEDURE MINFIT (INTEGER VALUE N; LONG REAL VALUE EPS, TOL;
+ LONG REAL ARRAY AB(*,*); LONG REAL ARRAY Q(*));
+ BEGIN COMMENT: AN IMPROVED VERSION OF MINFIT, SEE GOLUB &
+ REINSCH (1969), RESTRICTED TO M = N, P = 0.
+ THE SINGULAR VALUES OF THE ARRAY AB ARE
+ RETURNED IN Q, AND AB IS OVERWRITTEN WITH
+ THE ORTHOGONAL MATRIX V SUCH THAT
+ U.DIAG(Q) = AB.V,
+ WHERE U IS ANOTHER ORTHOGONAL MATRIX;
+ INTEGER L, KT;
+ LONG REAL C,F,G,H,S,X,Y,Z;
+ LONG REAL ARRAY E(1::N);
+ COMMENT: HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM;
+ G := X := 0;
+ FOR I := 1 UNTIL N DO
+ BEGIN
+ E(I) := G; S := 0; L := I+1;
+ FOR J := I UNTIL N DO S := S+AB(J,I)**2;
+ IF S<TOL THEN G := 0 ELSE
+ BEGIN
+ F := AB(I,I); G := IF F<0 THEN LONGSQRT(S)
+ ELSE -LONGSQRT(S);
+ H := F*G-S; AB(I,I) := F-G;
+ FOR J := L UNTIL N DO
+ BEGIN F := 0;
+ FOR K := I UNTIL N DO F := F + AB(K,I)*AB(K,J);
+ F := F/H;
+ FOR K := I UNTIL N DO AB(K,J) := AB(K,J) + F*AB(K,I)
+ END J
+ END S;
+ Q(I):=G; S:=0;
+ IF I<=N THEN FOR J := L UNTIL N DO
+ S:=S+AB(I,J)**2;
+ IF S<TOL THEN G := 0 ELSE
+ BEGIN
+ F := AB(I,I+1); G := IF F<0 THEN LONGSQRT(S)
+ ELSE -LONGSQRT(S);
+ H := F*G-S; AB(I,I+1) := F - G;
+ FOR J := L UNTIL N DO E(J) := AB(I,J)/H;
+ FOR J := L UNTIL N DO
+ BEGIN S := 0;
+ FOR K := L UNTIL N DO S := S + AB(J,K)*AB(I,K);
+ FOR K := L UNTIL N DO AB(J,K) := AB(J,K) + S*E(K)
+ END J
+ END S;
+ Y := ABS(Q(I)) + ABS(E(I)) ; IF Y >X THEN X := Y
+ END I;
+
+ COMMENT: ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS;
+ FOR I := N STEP -1 UNTIL 1 DO
+ BEGIN
+ IF G not =0 THEN
+ BEGIN
+ H := AB(I,I+1)*G;
+ FOR J := L UNTIL N DO AB(J,I) := AB(I,J)/H;
+ FOR J := L UNTIL N DO
+ BEGIN S := 0;
+ FOR K := L UNTIL N DO S := S + AB(I,K)*AB(K,J);
+ FOR K := L UNTIL N DO AB(K,J) := AB(K,J) + S*AB(K,I)
+ END J
+ END G;
+ FOR J := L UNTIL N DO AB(I,J) := AB(J,I) := 0;
+ AB(I,I) := 1; G := E(I); L := I
+ END I;
+
+ COMMENT: DIAGONALIZATION OF THE BIDIAGONAL FORM;
+ EPS := EPS*X;
+ FOR K := N STEP -1 UNTIL 1 DO
+ BEGIN KT := 0;
+ TESTFSPLITTING:
+ KT := KT + 1; IF KT > 30 THEN
+ BEGIN E(K) := 0L;
+ WRITE ("QR FAILED")
+ END;
+ FOR L2 := K STEP -1 UNTIL 1 DO
+ BEGIN
+ L := L2;
+ IF ABS(E(L))<=EPS THEN GOTO TESTFCONVERGENCE;
+ IF ABS(Q(L-1))<=EPS THEN GOTO CANCELLATION
+ END L2;
+
+ COMMENT: CANCELLATION OF E(L) IF L>1;
+ CANCELLATION:
+ C := 0; S := 1;
+ FOR I := L UNTIL K DO
+ BEGIN
+ F := S*E(I); E(I) := C*E(I);
+ IF ABS(F)<=EPS THEN GOTO TESTFCONVERGENCE;
+ G := Q(I); Q(I) := H := IF ABS(F) < ABS(G) THEN
+ ABS(G)*LONGSQRT(1 + (F/G)**2) ELSE IF F = 0 THEN
+ ABS(F)*LONGSQRT(1 + (G/F)**2) ELSE 0;
+ IF H = 0 THEN G := H := 1;
+ COMMENT: THE ABOVE REPLACES Q(I):=H:=LONGSQRT(G*G+F*F)
+ WHICH MAY GIVE INCORRECT RESULTS IF THE
+ SQUARES UNDERFLOW OR IF F = G = 0;
+ C := G/H; S := -F/H
+ END I;
+
+ TESTFCONVERGENCE:
+ Z := Q(K); IF L=K THEN GOTO CONVERGENCE;
+
+ COMMENT: SHIFT FROM BOTTOM 2*2 MINOR;
+ X := Q(L); Y := Q(K-1); G := E(K-1); H := E(K);
+ F := ((Y-Z)*(Y+Z) + (G-H)*(G+H))/(2*H*Y);
+ G := LONGSQRT(F*F+1);
+ F := ((X-Z)*(X+Z)+H*(Y/(IF F<0 THEN F-G ELSE F+G)-H))/X;
+
+ COMMENT: NEXT QR TRANSFORMATION;
+ C := S := 1;
+ FOR I := L+1 UNTIL K DO
+ BEGIN
+ G := E(I); Y := Q(I); H := S*G; G := G*C;
+ E(I-1) := Z := IF ABS(F) < ABS(H) THEN
+ ABS(H)*LONGSQRT(1 + (F/H)**2) ELSE IF F not = 0 THEN
+ ABS(F)*LONGSQRT(1 + (H/F)**2) ELSE 0;
+ IF Z = 0 THEN Z := F := 1 ;
+ C := F/Z; S := H/Z;
+ F := X*C + G*S; G := -X*S +G*C; H := Y*S;
+ Y := Y*C;
+ FOR J := 1 UNTIL N DO
+ BEGIN
+ X := AB(J,I-1); Z := AB(J,I);
+ AB(J,I-1) := X*C + Z*S; AB(J,I) := -X*S + Z*C
+ END J;
+ Q(I-1) := Z := IF ABS(F) < ABS(H) THEN ABS(H)*
+ LONGSQRT (1 + (F/H)**2) ELSE IF F not = 0 THEN
+ ABS(F)*LONGSQRT(1 + (H/F)**2) ELSE 0;
+ IF Z = 0 THEN Z := F := 1;
+ C := F/Z; S := H/Z;
+ F := C*G + S*Y; X := -S*G + C*Y
+ END I ;
+ E(L) := 0; E(K) := F; Q(K) := X;
+ GO TO TESTFSPLITTING;
+
+ CONVERGENCE:
+ IF Z<0 THEN
+ BEGIN COMMENT: Q(K) IS MADE NON-NEG;
+ Q(K) := -Z;
+ FOR J := 1 UNTIL N DO AB(J,K) := -AB(J,K)
+ END Z
+ END K
+ END MINFIT;
+
+ PROCEDURE SORT;
+ BEGIN COMMENT: SORTS THE ELEMENTS OF D AND CORRESPONDING
+ COLUMNS OF V INTO DESCENDING ORDER;
+ INTEGER K;
+ LONG REAL S;
+ FOR I := 1 UNTIL N - 1 DO
+ BEGIN K := I; S := D(I); FOR J := I + 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 UNTIL N DO
+ BEGIN S := V(J,I); V(J,I) := V(J,K); V(J,K) := S
+ END
+ END
+ END
+ END SORT;
+
+
+ PROCEDURE MATPRINT (STRING(80) VALUE S; LONG REAL ARRAY
+ V(*,*); INTEGER VALUE M, N);
+ BEGIN COMMENT: PRINTS M X N MATRIX V COLUMN BY COLUMN;
+ WRITE (S);
+ FOR K := 1 UNTIL (N + 7) DIV 8 DO
+ BEGIN FOR I := 1 UNTIL M DO
+ BEGIN IOCONTROL(2);
+ FOR J := 8*K - 7 UNTIL (IF N < (8*K) THEN N ELSE 8*K)
+ DO WRITEON (ROUNDTOREAL (V (I,J)))
+ END;
+ WRITE (" "); IOCONTROL(2)
+ END
+ END MATPRINT;
+
+ PROCEDURE VECPRINT (STRING(32) VALUE S; LONG REAL ARRAY V(*);
+ INTEGER VALUE N);
+ BEGIN COMMENT: PRINTS THE HEADING S AND N-VECTOR V;
+ WRITE(S);
+ FOR I := 1 UNTIL N DO WRITEON(ROUNDTOREAL(V(I)))
+ END VECPRINT;
+
+ PROCEDURE MIN (INTEGER VALUE J, NITS; LONG REAL VALUE
+ RESULT D2, X1; LONG REAL VALUE F1; LOGICAL VALUE FK);
+ BEGIN COMMENT:
+ MINIMIZES F FROM X IN THE DIRECTION V(*,J)
+ UNLESS J<1, WHEN A QUADRATIC SEARCH IS DONE
+ IN THE PLANE DEFINED BY Q0, Q1 AND X.
+ D2 AN APPROXIMATION TO HALF F'' (OR ZERO),
+ X1 AN ESTIMATE OF DISTANCE TO MINIMUM,
+ RETURNED AS THE DISTANCE FOUND.
+ IF FK = TRUE THEN F1 IS FLIN(X1), OTHERWISE
+ X1 AND F1 ARE IGNORED ON ENTRY UNLESS FINAL
+ FX > F1. NITS CONTROLS THE NUMBER OF TIMES
+ AN ATTEMPT IS MADE TO HALVE THE INTERVAL.
+ SIDE EFFECTS: USES AND ALTERS X, FX, NF, NL.
+ IF J < 1 USES VARIABLES Q... .
+ USES H, N, T, M2, M4, LDT, DMIN, MACHEPS;
+
+ LONG REAL PROCEDURE FLIN (LONG REAL VALUE L);
+ COMMENT: THE FUNCTION OF ONE VARIABLE L WHICH IS
+ MINIMIZED BY PROCEDURE MIN;
+ BEGIN LONG REAL ARRAY T(1::N);
+ IF J > 0 THEN
+ BEGIN COMMENT: LINEAR SEARCH;
+ FOR I := 1 UNTIL N DO T(I) := X(I) + L*V(I,J)
+ END
+ ELSE
+ BEGIN COMMENT: SEARCH ALONG A 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 UNTIL N DO T(I) := QA*Q0(I) + QB*X(I) + QC*Q1(I)
+ END;
+ COMMENT: INCREMENT FUNCTION EVALUATION COUNTER;
+ NF := NF + 1;
+ F(T,N)
+ END FLIN;
+
+ INTEGER K; LOGICAL DZ;
+ LONG REAL X2, XM, F0, F2, FM, D1, T2, S, SF1, SX1;
+ SF1 := F1; SX1 := X1;
+ K := 0; XM := 0; F0 := FM := FX; DZ := (D2 < MACHEPS);
+ COMMENT: FIND STEP SIZE;
+ S := 0; FOR I := 1 UNTIL N DO S := S + X(I)**2;
+ S := LONGSQRT(S);
+ T2:= M4*LONGSQRT(ABS(FX)/(IF DZ THEN DMIN ELSE D2)
+ + S*LDT) + M2*LDT;
+ S := M4*S + T;
+ 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 not FK OR (ABS(X1) < T2) THEN
+ BEGIN X1 := IF X1 >= 0L THEN T2 ELSE -T2;
+ F1 := FLIN(X1)
+ END;
+ IF F1 <= FM THEN BEGIN XM := X1; FM := F1 END;
+ L0: IF DZ THEN
+ BEGIN COMMENT: EVALUATE FLIN AT ANOTHER POINT AND
+ ESTIMATE THE SECONO DERIVATIVE;
+ X2 := IF F0 < F1 THEN -X1 ELSE 2*X1;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;
+ COMMENT: ESTIMATE FIRST DERIVATIVE AT 0;
+ D1 := (F1 - F0)/X1 - X1*D2; DZ := TRUE;
+ COMMENT: PREDICT MINIMUM;
+ X2 := IF D2 <- SMALL THEN (IF D1 < 0 THEN H ELSE -H) ELSE
+ -0.5L*D1/D2;
+ IF ABS(X2) > H THEN X2 := IF X2 > 0 THEN H ELSE -H;
+ COMMENT: EVALUATE F AT THE PREDICTED M(NIMUM;
+ L1: F2 := FLIN(X2);
+ IF (K < NITS) AND (F2 > F0) THEN
+ BEGIN COMMENT: NO SUCCESS SO TRY AGAIN; K := K + 1;
+ IF (F0 < F1) AND ((X1*X2) > 0) THEN GO TO L0;
+ X2 := 0.5L*X2; GO TO L1
+ END;
+ COMMENT: INCREMENT ONE-DIMENSIONAL SEARCH COUNTER;
+ NL := NL + 1;
+ IF F2 > FM THEN X2 := XM ELSE FM := F2;
+ COMMENT: GET NEW ESTIMATE OF SECUND DERIVATIVE;
+ D2 := IF ABS(X2*(X2 - X1)) > SMALL THEN
+ (X2*(F1 - F0) - X1*(FM - F0))/(X1*X2*(X1 - X2))
+ ELSE IF K > 0 THEN 0 ELSE D2;
+ IF D2 <= SMALL THEN D2 := SMALL;
+ X1 := X2; FX := FM;
+ IF SF1 < FX THEN BEGIN FX := SF1; X1 := SX1 END;
+ COMMENT: UPDATE X FOR LINEAR SEARCH BUT NOT FOR PARABOLIC
+ PARABOLIC SEARCH;
+ IF J > 0 THEN FOR I := 1 UNTIL N DO X(I) := X(I) + X1*V(I,J)
+ END MIN;
+
+ PROCEDURE QUAD;
+ BEGIN COMMENT: LOOKS FOR THE MINIMUM ALONG A CURVE
+ DEFINED BY Q0, Q1 AND X;
+ LONG REAL L, S;
+ S := FX; FX := QF1; QF1 := S; QD1 := 0;
+ FOR I := 1 UNTIL N DO
+ BEGIN S := X(I); X(I) := L := Q1(I); Q1(I):= S;
+ QD1 := QD1 + (S - L)**2
+ END;
+ L := QD1 := LONGSQRT(QD1); S := 0;
+ IF (QD0 > 0) AND (QD1 > 0) 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 := 1 UNTIL N DO
+ BEGIN S := Q0(I); Q0(1) := X(I);
+ X(I) := QA*S + QB*X(I) + QC*Q1(I)
+ END
+ END QUAD;
+
+ PROCEDURE PRINT;
+ COMMENT: THE VARIABLE FMIN IS GLOBAL, AND ESTIMATES THE
+ VALUE OF F AT THE MINIMUM: USED ONLY FOR
+ PRINTING LOG(FX - FMIN);
+ IF PRIN > 0 THEN
+ BEGIN INTEGER SVINT; long real fmin;
+ SVINT := I_W;
+ I_W := 10; % print integers in 10 column fields %
+ WRITE (NL, NF, FX);
+ COMMENT: IF THE NEXT TWO LINES ARE OMITTED THEN FMIN IS
+ NOT REQUIRED;
+ IF FX <= FMIN THEN WRITEON (" UNDEFINED ") ELSE
+ WRITEON (ROUNDTOREAL (LONGLOG (FX - FMIN )));
+ COMMENT: "IOCONTROL(2)" MOVES TO THE NEXT LINE;
+ IF N > 4 THEN IOCONTROL(2);
+ IF (N <= 4) OR (PRIN > 2) THEN
+ FOR I := 1 UNTIL N DO WRITEON(ROUNDTOREAL(X(I)));
+ IOCONTROL(2);
+ I_W := SVINT
+ END PRINT;
+
+ LOGICAL ILLC;
+ INTEGER NL, NF, KL, KT, KTM;
+ LONG 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;
+ LONG REAL ARRAY D, Y, Z, Q0, Q1 (1::N);
+ LONG REAL ARRAY V (1::N, 1::N);
+
+ COMMENT: INITIALIZATION;
+ COMMENT: MACHINE DEPENDENT NUMBERS;
+ SMALL := MACHEPS**2; VSMALL := SMALL**2;
+ LARGE := 1L/SMALL; VLARGE := 1L/VSMALL;
+ M2 := LONGSQRT(MACHEPS); M4 := LONGSQRT(M2);
+
+ COMMENT: HEURISTIC NUMBERS
+ •••••••••••••
+
+ IF AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
+ POSSIBLE! THEN SET SCBD := 10, OTHERWISE 1,
+ IF THE PROBLEM IS KNOWN TO BE ILLCONDITIONED SET
+ ILLC := TRUE, OTHERWISE FALSE,
+ KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE
+ THE ALGORITHM TERMINATES (SEE SECTION 6). KTM = 4, IS VERY
+ CAUTIOUS: USUALLY KTM = 1 IS SATISFACTORY;
+
+ SCBD := 1; ILLC := FALSE; KTM := 1;
+
+ LDFAC := IF ILLC THEN 0.1 ELSE 0.01;
+ KT := NL := 0; NF := 1; QF1 := FX := F(X,N);
+ T := T2 := SMALL + ABS(T); DMIN := SMALL;
+ IF H < (100*T) THEN H := 100*T; LDT := H;
+ FOR I := 1 UNTIL N DO FOR J := 1 UNTIL N DO
+ V(I,J) := IF I = J THEN 1L ELSE 0L;
+ D(1) := QD0 := 0; FOR I := 1 UNTIL N DO Q1(I) := X(I);
+ PRINT;
+
+ COMMENT: MAIN LOOP;
+ L0: SF := D(1); D(1) := S := 0;
+ COMMENT: MINIMIZE ALONG FIRST DIRECTION;
+ MIN (1, 2, D(1), S, FX, FALSE);
+ IF S <= 0 THEN FOR I := 1 UNTIL N DO V(I,1) := -V(I,1);
+ IF (SF <= (0.9*D(1))) OR ((0.9*SF) >= D(1)) THEN
+ FOR I := 2 UNTIL N DO D(I) := 0;
+ FOR K := 2 UNTIL N DO
+ BEGIN FOR I := 1 UNTIL N DO Y(I) := X(I); SF := FX;
+ ILLC := ILLC OR (KT > 0);
+ L1: KL := K; DF := 0; IF ILLC THEN
+ BEGIN COMMENT: RANDOM STEP TO GET OFF RESOLUTION VALLEY;
+ FOR I := 1 UNTIL N DO
+ BEGIN S := Z(I) := (0.1*LDT + T2*10**KT)*(RANDOM(I)-0.5L);
+ COMMENT: PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM
+ NUMBER UNIFORMLY DISTRIBUTED IN (0, 1) AND
+ THAT ANY INITIALIZATION OF THE RANDOM NUMBER
+ GENERATOR HAS ALREADY BEEN DONE;
+ FOR J := 1 UNTIL N DO X(J) := X(J) + S*V(J,I)
+ END;
+ FX := F(X,N); NF := NF + 1
+ END;
+ FOR K2 := K UNTIL N DO
+ BEGIN SL := FX; S := 0;
+ COMMENT: MINIMIZE ALONG "NON-CONJUGATE" DIRECTIONS;
+ 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 not ILLC AND (DF < ABS( 100*MACHEPS*FX)) THEN
+ BEGIN COMMENT: NO SUCCESS ILLC = FALSE SO TRY ONCE
+ WITH ILLC = TRUE;
+ ILLC := TRUE; GO TO L1
+ END;
+ IF (K = 2) AND (PRIN > 1) THEN VECPRINT ("NEW D", D, N);
+ FOR K2 := 1 UNTIL K - 1 DO
+ BEGIN COMMENT: MINIMIZE ALONG "CONJUGATE" DIRECTIONS;
+ S := 0; MIN (K2, 2, D(K2), S, FX, FALSE)
+ END;
+ F1 := FX; FX := SF; LDS := 0;
+ FOR I := 1 UNTIL N DO
+ BEGIN SL := X(I); X(I) := Y(I); SL := Y(I) := SL - Y(I);
+ LDS := LDS + SL*SL
+ END;
+ LDS := LONGSQRT(LDS); IF LDS > SMALL THEN
+ BEGIN COMMENT: THROW AWAY DIRECTION KL AND MINIMIZE
+ ALONG THE NEW "CONJUGATE" DIRECTION;
+ FOR I := KL - 1 STEP -1 UNTIL K DO
+ BEGIN FOR J := 1 UNTIL N DO V(J,I + 1) := V(J,I);
+ D(I + 1) := D(I)
+ END;
+ D(K) := 0; FOR I := 1 UNTIL N DO V(I,K) := Y(I)/LDS;
+ MIN (K, 4, D(K), LDS, F1, TRUE);
+ IF LDS <= 0 THEN
+ BEGIN LDS := -LDS;
+ FOR I := 1 UNTIL N DO V(I,K) := -V(I,K)
+ END
+ END;
+ LDT := LDFAC*LDT; IF LDT < LDS THEN LDT := LDS;
+ PRINT;
+ T2 := 0; FOR I := 1 UNTIL N DO T2 := T2 + X(I)**2;
+ T2 := M2*LONGSQRT(T2) + T;
+ COMMENT: SEE IF STEP LENGTH EXCEEDS HALF THE TOLERANCE;
+ KT := IF LDT > (0.5*T2) THEN 0 ELSE KT + 1;
+ IF KT > KTM THEN GO TO L2
+ END;
+ COMMENT: TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE STUCK
+ IN A CURVED VALLEY;
+ QUAD;
+ DN := 0; FOR I := 1 UNTIL N DO
+ BEGIN D(I) := 1/LONGSQRT(D(I));
+ IF DN < D(I) THEN DN : = D(I)
+ END;
+ IF PRIN > 3 THEN MATPRINT ("NEW DIRECTIONS", V, N, N);
+ FOR J := 1 UNTIL N DO
+ BEGIN S := D(J)/DN;
+ FOR I := 1 UNTIL N DO V(I,J) := S*V(I,J)
+ END;
+ IF SCBD > 1 THEN
+ BEGIN COMMENT: SCALE AXES TO TRY TO REDUCE CONDITION
+ NUMBER;
+ S := VLARGE; FOR I := 1 UNTIL N DO
+ BEGIN SL := 0; FOR J := 1 UNTIL N DO SL := SL+V(I,J)**2;
+ Z(I) := LONGSQRT(SL);
+ IF Z(I) < M4 THEN Z(I) := M4; IF S > Z(I) THEN S := Z(I)
+ END;
+ FOR I := 1 UNTIL N DO
+ BEGIN SL := S/Z(I); Z(I) := 1/SL; IF Z(I) > SCBD THEN
+ BEGIN SL := 1/SCBD; Z(I) := SCBD
+ END;
+ FOR J := 1 UNTIL N DO V(I,J) := SL*V(I,J)
+ END
+ END;
+ COMMENT: TRANSPOSE V FOR MINFIT LINE BEFORE WAS OMMITTED IN PUBLICATION;
+ FOR I := 2 UNTIL N DO FOR J := 1 UNTIL I - 1 DO
+ BEGIN S := V(I,J); V(I,J) := V(J,I); V(J,I) := S END;
+ COMMENT: FIND THE SINGULAR VALUE DECOMPOSITION OF V, THIS
+ GIVES THE EIGENVALUES AND PRINCIPAL AXES OF THE
+ APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE
+ CONDITION NUMBER;
+ MINFIT (N, MACHEPS, VSMALL, V, D);
+ IF SCBD > 1 THEN
+ BEGIN COMMENT: UNSCALlNG; FOR I := 1 UNTIL N DO
+ BEGIN S := Z(I) ;
+ FOR J := 1 UNTIL N DO V(I,J) := S*V(I,J)
+ END;
+ FOR I := 1 UNTIL N DO
+ BEGIN S := 0; FOR J := 1 UNTIL N DO S := S + V(J,I)**2;
+ S := LONGSQRT(S); D(I) := S*D(I); S := 1/S;
+ FOR J := 1 UNTIL N DO V(J,I) := S*V(J,I)
+ END
+ END;
+ FOR I := 1 UNTIL N DO
+ BEGIN D(I) := IF (DN*D(I)) > LARGE THEN VSMALL ELSE
+ IF (DN*D(I)) < SMALL THEN VLARGE ELSE (DN*D(I))**(-2)
+ END;
+ COMMENT: SORT NEW EIGENVALUES AND EIGENVECTORS;
+ SORT;
+ DMIN := D(N) ; IF DMIN < SMALL THEN DMIN := SMALL;
+ ILLC := (M2*D(1)) > DMIN;
+ IF (PRIN > 1) AND (SCBD > 1) THEN
+ VECPRINT ("SCALE FACTORS", Z, N);
+ IF PRIN > 1 THEN VECPRINT ("EIGENVALUES OF A", D, N);
+ IF PRIN > 3 THEN MATPRINT ("EIGENVECTORS OF A", V, N, N);
+ COMMENT: GO BACK TO MAIN LOOP;
+ GO TO L0;
+ L2: IF PRIN > 0 THEN VECPRINT ("X IS", X, N);
+ FX
+ END ALGOLPRAXIS.
git://henry.ined.fr / .git / commitdiff