From 1609c72bcd97ccf2707cec75b457792a53c1509a Mon Sep 17 00:00:00 2001 From: "N. Brouard" Date: Sat, 3 Feb 2024 23:35:51 +0000 Subject: [PATCH] Summary: praxis procedure as published by brent 1973 With the minimum of modifications to run the tests and to be called by a C main. --- src/praxis.alw | 523 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 523 insertions(+) create mode 100644 src/praxis.alw diff --git a/src/praxis.alw b/src/praxis.alw new file mode 100644 index 0000000..a849dca --- /dev/null +++ b/src/praxis.alw @@ -0,0 +1,523 @@ + 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 SX 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. -- 2.43.0