1: LONG REAL PROCEDURE ALGOLPRAXIS (LONG REAL VALUE T, MACHEPS, H ;
2: INTEGER VALUE N, PRIN;
3: LONG REAL ARRAY X(*);
4: LONG REAL PROCEDURE F(long real array x(*);
5: integer value n));
6: BEGIN COMMENT:
7: THIS PROCEDURE MINIMIZES THE FONCTION F(X, N) OF N
8: VARIABLES X(1), ... X(N), USING THE PRINCIPAL AXIS METHOD.
9: ON ENTRY X HOLDS A GUESS, ON RETURN IT HOLDS THE ESTIMATED
10: POINT OF MINIMUM, WITH (HOPEFULLY) |ERROR| <
11: SQRT(MACHEPS)*|X| + T, WHERE MACHEPS IS THE MACHINE
12: PRECISION, THE SMALLEST NUMBER SUCH THAT 1 + MACHEPS > 1,
13: T IS A TOLERANCE, AND |.| IS THE 2-NORM. H IS THE MAXIMUM
14: STEP SIZE: SET TO ABOUT THE MAXIMUM EXPECTED DISTANCE FROM
15: THE GUESS TO THE MINIMUM (IF H IS SET TOO SMALL OR TOO
16: LARGE THEN THE INITIAL RATE OF CONVERGENCE WILL BE SLOW).
17: THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS
18: AFTER PROCEDURE QUAD.
19: PRIN CONTROLS THE PRINTING OF INTERMEDIATE RESULTS.
20: IF PRIN = 0, NO RESULTS ARE PRINTED.
21: IF PRIN = 1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR
22: MINIMIZATIONS, AND FINAL X IS PRINTED, BUT INTERMEDIATE
23: X ONLY IF N <= 4.
24: IF PRIN = 2, EIGENVALUES OF A AND SCALE FACTORS ARE ALSO PRINTED.
25: IF PRIN = 3, F AND X ARE PRINTED AFTER EVERY FEW LINEAR MINIMIZATIONS.
26: IF PRIN = 4, EIGENVECTORS ARE ALSO PRINTED.
27: FMIN IS A GLOBAL VARIABLE: SEE PROCEDURE PRINT.
28: RANDOM IS A PARAMETERLESS LONG REAL PROCEDURE WHICH RETURNS
29: A RANDOM NUMBER UNIFORMLY DISTRIBUTED IN (0, 1). ANY
30: INITIALIZATION MUST BE DONE BEFORE THE CALL TO PRAXIS.
31: THE PROCEDURE IS MACHINE-INDEPENDENT, APART FROM THE OUTPUT
32: STATEMENTS AND THE SPECIFICATION OF MACHEPS. WE ASSUME THAT
33: MACHEPS**(—4) DOES NOT OVERFLOW (IF IT DOES THEN MACHEPS MUST
34: BE INCREASED), AND THAT ON FLOATING-POINT UNDERFLOW THE
35: RESULT IS SET TO ZERO;
36:
37: LONG REAL PROCEDURE RANDOM(INTEGER VALUE NAUGHT);
38: ALGOL "random";
39:
40: PROCEDURE MINFIT (INTEGER VALUE N; LONG REAL VALUE EPS, TOL;
41: LONG REAL ARRAY AB(*,*); LONG REAL ARRAY Q(*));
42: BEGIN COMMENT: AN IMPROVED VERSION OF MINFIT, SEE GOLUB &
43: REINSCH (1969), RESTRICTED TO M = N, P = 0.
44: THE SINGULAR VALUES OF THE ARRAY AB ARE
45: RETURNED IN Q, AND AB IS OVERWRITTEN WITH
46: THE ORTHOGONAL MATRIX V SUCH THAT
47: U.DIAG(Q) = AB.V,
48: WHERE U IS ANOTHER ORTHOGONAL MATRIX;
49: INTEGER L, KT;
50: LONG REAL C,F,G,H,S,X,Y,Z;
51: LONG REAL ARRAY E(1::N);
52: COMMENT: HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM;
53: G := X := 0;
54: FOR I := 1 UNTIL N DO
55: BEGIN
56: E(I) := G; S := 0; L := I+1;
57: FOR J := I UNTIL N DO S := S+AB(J,I)**2;
58: IF S<TOL THEN G := 0 ELSE
59: BEGIN
60: F := AB(I,I); G := IF F<0 THEN LONGSQRT(S)
61: ELSE -LONGSQRT(S);
62: H := F*G-S; AB(I,I) := F-G;
63: FOR J := L UNTIL N DO
64: BEGIN F := 0;
65: FOR K := I UNTIL N DO F := F + AB(K,I)*AB(K,J);
66: F := F/H;
67: FOR K := I UNTIL N DO AB(K,J) := AB(K,J) + F*AB(K,I)
68: END J
69: END S;
70: Q(I):=G; S:=0;
71: IF I<=N THEN FOR J := L UNTIL N DO
72: S:=S+AB(I,J)**2;
73: IF S<TOL THEN G := 0 ELSE
74: BEGIN
75: F := AB(I,I+1); G := IF F<0 THEN LONGSQRT(S)
76: ELSE -LONGSQRT(S);
77: H := F*G-S; AB(I,I+1) := F - G;
78: FOR J := L UNTIL N DO E(J) := AB(I,J)/H;
79: FOR J := L UNTIL N DO
80: BEGIN S := 0;
81: FOR K := L UNTIL N DO S := S + AB(J,K)*AB(I,K);
82: FOR K := L UNTIL N DO AB(J,K) := AB(J,K) + S*E(K)
83: END J
84: END S;
85: Y := ABS(Q(I)) + ABS(E(I)) ; IF Y >X THEN X := Y
86: END I;
87:
88: COMMENT: ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS;
89: FOR I := N STEP -1 UNTIL 1 DO
90: BEGIN
91: IF G not =0 THEN
92: BEGIN
93: H := AB(I,I+1)*G;
94: FOR J := L UNTIL N DO AB(J,I) := AB(I,J)/H;
95: FOR J := L UNTIL N DO
96: BEGIN S := 0;
97: FOR K := L UNTIL N DO S := S + AB(I,K)*AB(K,J);
98: FOR K := L UNTIL N DO AB(K,J) := AB(K,J) + S*AB(K,I)
99: END J
100: END G;
101: FOR J := L UNTIL N DO AB(I,J) := AB(J,I) := 0;
102: AB(I,I) := 1; G := E(I); L := I
103: END I;
104:
105: COMMENT: DIAGONALIZATION OF THE BIDIAGONAL FORM;
106: EPS := EPS*X;
107: FOR K := N STEP -1 UNTIL 1 DO
108: BEGIN KT := 0;
109: TESTFSPLITTING:
110: KT := KT + 1; IF KT > 30 THEN
111: BEGIN E(K) := 0L;
112: WRITE ("QR FAILED")
113: END;
114: FOR L2 := K STEP -1 UNTIL 1 DO
115: BEGIN
116: L := L2;
117: IF ABS(E(L))<=EPS THEN GOTO TESTFCONVERGENCE;
118: IF ABS(Q(L-1))<=EPS THEN GOTO CANCELLATION
119: END L2;
120:
121: COMMENT: CANCELLATION OF E(L) IF L>1;
122: CANCELLATION:
123: C := 0; S := 1;
124: FOR I := L UNTIL K DO
125: BEGIN
126: F := S*E(I); E(I) := C*E(I);
127: IF ABS(F)<=EPS THEN GOTO TESTFCONVERGENCE;
128: G := Q(I); Q(I) := H := IF ABS(F) < ABS(G) THEN
129: ABS(G)*LONGSQRT(1 + (F/G)**2) ELSE IF F = 0 THEN
130: ABS(F)*LONGSQRT(1 + (G/F)**2) ELSE 0;
131: IF H = 0 THEN G := H := 1;
132: COMMENT: THE ABOVE REPLACES Q(I):=H:=LONGSQRT(G*G+F*F)
133: WHICH MAY GIVE INCORRECT RESULTS IF THE
134: SQUARES UNDERFLOW OR IF F = G = 0;
135: C := G/H; S := -F/H
136: END I;
137:
138: TESTFCONVERGENCE:
139: Z := Q(K); IF L=K THEN GOTO CONVERGENCE;
140:
141: COMMENT: SHIFT FROM BOTTOM 2*2 MINOR;
142: X := Q(L); Y := Q(K-1); G := E(K-1); H := E(K);
143: F := ((Y-Z)*(Y+Z) + (G-H)*(G+H))/(2*H*Y);
144: G := LONGSQRT(F*F+1);
145: F := ((X-Z)*(X+Z)+H*(Y/(IF F<0 THEN F-G ELSE F+G)-H))/X;
146:
147: COMMENT: NEXT QR TRANSFORMATION;
148: C := S := 1;
149: FOR I := L+1 UNTIL K DO
150: BEGIN
151: G := E(I); Y := Q(I); H := S*G; G := G*C;
152: E(I-1) := Z := IF ABS(F) < ABS(H) THEN
153: ABS(H)*LONGSQRT(1 + (F/H)**2) ELSE IF F not = 0 THEN
154: ABS(F)*LONGSQRT(1 + (H/F)**2) ELSE 0;
155: IF Z = 0 THEN Z := F := 1 ;
156: C := F/Z; S := H/Z;
157: F := X*C + G*S; G := -X*S +G*C; H := Y*S;
158: Y := Y*C;
159: FOR J := 1 UNTIL N DO
160: BEGIN
161: X := AB(J,I-1); Z := AB(J,I);
162: AB(J,I-1) := X*C + Z*S; AB(J,I) := -X*S + Z*C
163: END J;
164: Q(I-1) := Z := IF ABS(F) < ABS(H) THEN ABS(H)*
165: LONGSQRT (1 + (F/H)**2) ELSE IF F not = 0 THEN
166: ABS(F)*LONGSQRT(1 + (H/F)**2) ELSE 0;
167: IF Z = 0 THEN Z := F := 1;
168: C := F/Z; S := H/Z;
169: F := C*G + S*Y; X := -S*G + C*Y
170: END I ;
171: E(L) := 0; E(K) := F; Q(K) := X;
172: GO TO TESTFSPLITTING;
173:
174: CONVERGENCE:
175: IF Z<0 THEN
176: BEGIN COMMENT: Q(K) IS MADE NON-NEG;
177: Q(K) := -Z;
178: FOR J := 1 UNTIL N DO AB(J,K) := -AB(J,K)
179: END Z
180: END K
181: END MINFIT;
182:
183: PROCEDURE SORT;
184: BEGIN COMMENT: SORTS THE ELEMENTS OF D AND CORRESPONDING
185: COLUMNS OF V INTO DESCENDING ORDER;
186: INTEGER K;
187: LONG REAL S;
188: FOR I := 1 UNTIL N - 1 DO
189: BEGIN K := I; S := D(I); FOR J := I + 1 UNTIL N DO
190: IF D(J) > S THEN
191: BEGIN K := J; S := D(J) END;
192: IF K > I THEN
193: BEGIN D(K) := D(I); D(I) := S; FOR J := 1 UNTIL N DO
194: BEGIN S := V(J,I); V(J,I) := V(J,K); V(J,K) := S
195: END
196: END
197: END
198: END SORT;
199:
200:
201: PROCEDURE MATPRINT (STRING(80) VALUE S; LONG REAL ARRAY
202: V(*,*); INTEGER VALUE M, N);
203: BEGIN COMMENT: PRINTS M X N MATRIX V COLUMN BY COLUMN;
204: WRITE (S);
205: FOR K := 1 UNTIL (N + 7) DIV 8 DO
206: BEGIN FOR I := 1 UNTIL M DO
207: BEGIN IOCONTROL(2);
208: FOR J := 8*K - 7 UNTIL (IF N < (8*K) THEN N ELSE 8*K)
209: DO WRITEON (ROUNDTOREAL (V (I,J)))
210: END;
211: WRITE (" "); IOCONTROL(2)
212: END
213: END MATPRINT;
214:
215: PROCEDURE VECPRINT (STRING(32) VALUE S; LONG REAL ARRAY V(*);
216: INTEGER VALUE N);
217: BEGIN COMMENT: PRINTS THE HEADING S AND N-VECTOR V;
218: WRITE(S);
219: FOR I := 1 UNTIL N DO WRITEON(ROUNDTOREAL(V(I)))
220: END VECPRINT;
221:
222: PROCEDURE MIN (INTEGER VALUE J, NITS; LONG REAL VALUE
223: RESULT D2, X1; LONG REAL VALUE F1; LOGICAL VALUE FK);
224: BEGIN COMMENT:
225: MINIMIZES F FROM X IN THE DIRECTION V(*,J)
226: UNLESS J<1, WHEN A QUADRATIC SEARCH IS DONE
227: IN THE PLANE DEFINED BY Q0, Q1 AND X.
228: D2 AN APPROXIMATION TO HALF F'' (OR ZERO),
229: X1 AN ESTIMATE OF DISTANCE TO MINIMUM,
230: RETURNED AS THE DISTANCE FOUND.
231: IF FK = TRUE THEN F1 IS FLIN(X1), OTHERWISE
232: X1 AND F1 ARE IGNORED ON ENTRY UNLESS FINAL
233: FX > F1. NITS CONTROLS THE NUMBER OF TIMES
234: AN ATTEMPT IS MADE TO HALVE THE INTERVAL.
235: SIDE EFFECTS: USES AND ALTERS X, FX, NF, NL.
236: IF J < 1 USES VARIABLES Q... .
237: USES H, N, T, M2, M4, LDT, DMIN, MACHEPS;
238:
239: LONG REAL PROCEDURE FLIN (LONG REAL VALUE L);
240: COMMENT: THE FUNCTION OF ONE VARIABLE L WHICH IS
241: MINIMIZED BY PROCEDURE MIN;
242: BEGIN LONG REAL ARRAY T(1::N);
243: IF J > 0 THEN
244: BEGIN COMMENT: LINEAR SEARCH;
245: FOR I := 1 UNTIL N DO T(I) := X(I) + L*V(I,J)
246: END
247: ELSE
248: BEGIN COMMENT: SEARCH ALONG A PARABOLIC SPACE-CURVE;
249: QA := L*(L - QD1)/(QD0*(QD0 + QD1));
250: QB := (L + QD0)*(QD1 - L)/(QD0*QD1);
251: QC := L*(L + QD0)/(QD1*(QD0 + QD1));
252: FOR I := 1 UNTIL N DO T(I) := QA*Q0(I) + QB*X(I) + QC*Q1(I)
253: END;
254: COMMENT: INCREMENT FUNCTION EVALUATION COUNTER;
255: NF := NF + 1;
256: F(T,N)
257: END FLIN;
258:
259: INTEGER K; LOGICAL DZ;
260: LONG REAL X2, XM, F0, F2, FM, D1, T2, S, SF1, SX1;
261: SF1 := F1; SX1 := X1;
262: K := 0; XM := 0; F0 := FM := FX; DZ := (D2 < MACHEPS);
263: COMMENT: FIND STEP SIZE;
264: S := 0; FOR I := 1 UNTIL N DO S := S + X(I)**2;
265: S := LONGSQRT(S);
266: T2:= M4*LONGSQRT(ABS(FX)/(IF DZ THEN DMIN ELSE D2)
267: + S*LDT) + M2*LDT;
268: S := M4*S + T;
269: IF DZ AND (T2 > S) THEN T2 := S;
270: IF T2 < SMALL THEN T2 := SMALL;
271: IF T2 > (0.01*H) THEN T2 := 0.01*H;
272: IF FK AND (F1 <= FM) THEN BEGIN XM := X1; FM:=F1 END;
273: IF not FK OR (ABS(X1) < T2) THEN
274: BEGIN X1 := IF X1 >= 0L THEN T2 ELSE -T2;
275: F1 := FLIN(X1)
276: END;
277: IF F1 <= FM THEN BEGIN XM := X1; FM := F1 END;
278: L0: IF DZ THEN
279: BEGIN COMMENT: EVALUATE FLIN AT ANOTHER POINT AND
280: ESTIMATE THE SECONO DERIVATIVE;
281: X2 := IF F0 < F1 THEN -X1 ELSE 2*X1;F2:=FLIN(X2);
282: IF F2 <= FM THEN BEGIN XM := X2; FM := F2 END;
283: D2 := (X2*(F1 - F0) - X1*(F2 - F0))/(X1*X2*(X1 - X2))
284: END;
285: COMMENT: ESTIMATE FIRST DERIVATIVE AT 0;
286: D1 := (F1 - F0)/X1 - X1*D2; DZ := TRUE;
287: COMMENT: PREDICT MINIMUM;
288: X2 := IF D2 <- SMALL THEN (IF D1 < 0 THEN H ELSE -H) ELSE
289: -0.5L*D1/D2;
290: IF ABS(X2) > H THEN X2 := IF X2 > 0 THEN H ELSE -H;
291: COMMENT: EVALUATE F AT THE PREDICTED M(NIMUM;
292: L1: F2 := FLIN(X2);
293: IF (K < NITS) AND (F2 > F0) THEN
294: BEGIN COMMENT: NO SUCCESS SO TRY AGAIN; K := K + 1;
295: IF (F0 < F1) AND ((X1*X2) > 0) THEN GO TO L0;
296: X2 := 0.5L*X2; GO TO L1
297: END;
298: COMMENT: INCREMENT ONE-DIMENSIONAL SEARCH COUNTER;
299: NL := NL + 1;
300: IF F2 > FM THEN X2 := XM ELSE FM := F2;
301: COMMENT: GET NEW ESTIMATE OF SECUND DERIVATIVE;
302: D2 := IF ABS(X2*(X2 - X1)) > SMALL THEN
303: (X2*(F1 - F0) - X1*(FM - F0))/(X1*X2*(X1 - X2))
304: ELSE IF K > 0 THEN 0 ELSE D2;
305: IF D2 <= SMALL THEN D2 := SMALL;
306: X1 := X2; FX := FM;
307: IF SF1 < FX THEN BEGIN FX := SF1; X1 := SX1 END;
308: COMMENT: UPDATE X FOR LINEAR SEARCH BUT NOT FOR PARABOLIC
309: PARABOLIC SEARCH;
310: IF J > 0 THEN FOR I := 1 UNTIL N DO X(I) := X(I) + X1*V(I,J)
311: END MIN;
312:
313: PROCEDURE QUAD;
314: BEGIN COMMENT: LOOKS FOR THE MINIMUM ALONG A CURVE
315: DEFINED BY Q0, Q1 AND X;
316: LONG REAL L, S;
317: S := FX; FX := QF1; QF1 := S; QD1 := 0;
318: FOR I := 1 UNTIL N DO
319: BEGIN S := X(I); X(I) := L := Q1(I); Q1(I):= S;
320: QD1 := QD1 + (S - L)**2
321: END;
322: L := QD1 := LONGSQRT(QD1); S := 0;
323: IF (QD0 > 0) AND (QD1 > 0) AND (NL >= (3*N*N)) THEN
324: BEGIN MIN (0, 2, S, L, QF1, TRUE);
325: QA := L*(L - QD1)/(QD0*(QD0 + QD1));
326: QB := (L + QD0)*(QD1 - L)/(QD0*QD1);
327: QC := L*(L + QD0)/(QD1*(QD0 + QD1))
328: END
329: ELSE BEGIN FX := QF1; QA := QB := 0; QC := 1 END;
330: QD0 := QD1; FOR I := 1 UNTIL N DO
331: BEGIN S := Q0(I); Q0(1) := X(I);
332: X(I) := QA*S + QB*X(I) + QC*Q1(I)
333: END
334: END QUAD;
335:
336: PROCEDURE PRINT;
337: COMMENT: THE VARIABLE FMIN IS GLOBAL, AND ESTIMATES THE
338: VALUE OF F AT THE MINIMUM: USED ONLY FOR
339: PRINTING LOG(FX - FMIN);
340: IF PRIN > 0 THEN
341: BEGIN INTEGER SVINT; long real fmin;
342: SVINT := I_W;
343: I_W := 10; % print integers in 10 column fields %
344: WRITE (NL, NF, FX);
345: COMMENT: IF THE NEXT TWO LINES ARE OMITTED THEN FMIN IS
346: NOT REQUIRED;
347: IF FX <= FMIN THEN WRITEON (" UNDEFINED ") ELSE
348: WRITEON (ROUNDTOREAL (LONGLOG (FX - FMIN )));
349: COMMENT: "IOCONTROL(2)" MOVES TO THE NEXT LINE;
350: IF N > 4 THEN IOCONTROL(2);
351: IF (N <= 4) OR (PRIN > 2) THEN
352: FOR I := 1 UNTIL N DO WRITEON(ROUNDTOREAL(X(I)));
353: IOCONTROL(2);
354: I_W := SVINT
355: END PRINT;
356:
357: LOGICAL ILLC;
358: INTEGER NL, NF, KL, KT, KTM;
359: LONG REAL S, SL, DN, DMIN, FX, F1, LDS, LDT, SF, DF,
360: QF1, QD0, QD1, QA, QB, QC,
361: M2, M4, SMALL, VSMALL, LARGE, VLARGE, SCBD, LDFAC, T2;
362: LONG REAL ARRAY D, Y, Z, Q0, Q1 (1::N);
363: LONG REAL ARRAY V (1::N, 1::N);
364:
365: COMMENT: INITIALIZATION;
366: COMMENT: MACHINE DEPENDENT NUMBERS;
367: SMALL := MACHEPS**2; VSMALL := SMALL**2;
368: LARGE := 1L/SMALL; VLARGE := 1L/VSMALL;
369: M2 := LONGSQRT(MACHEPS); M4 := LONGSQRT(M2);
370:
371: COMMENT: HEURISTIC NUMBERS
372: •••••••••••••
373:
374: IF AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF
375: POSSIBLE! THEN SET SCBD := 10, OTHERWISE 1,
376: IF THE PROBLEM IS KNOWN TO BE ILLCONDITIONED SET
377: ILLC := TRUE, OTHERWISE FALSE,
378: KTM+1 IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE
379: THE ALGORITHM TERMINATES (SEE SECTION 6). KTM = 4, IS VERY
380: CAUTIOUS: USUALLY KTM = 1 IS SATISFACTORY;
381:
382: SCBD := 1; ILLC := FALSE; KTM := 1;
383:
384: LDFAC := IF ILLC THEN 0.1 ELSE 0.01;
385: KT := NL := 0; NF := 1; QF1 := FX := F(X,N);
386: T := T2 := SMALL + ABS(T); DMIN := SMALL;
387: IF H < (100*T) THEN H := 100*T; LDT := H;
388: FOR I := 1 UNTIL N DO FOR J := 1 UNTIL N DO
389: V(I,J) := IF I = J THEN 1L ELSE 0L;
390: D(1) := QD0 := 0; FOR I := 1 UNTIL N DO Q1(I) := X(I);
391: PRINT;
392:
393: COMMENT: MAIN LOOP;
394: L0: SF := D(1); D(1) := S := 0;
395: COMMENT: MINIMIZE ALONG FIRST DIRECTION;
396: MIN (1, 2, D(1), S, FX, FALSE);
397: IF S <= 0 THEN FOR I := 1 UNTIL N DO V(I,1) := -V(I,1);
398: IF (SF <= (0.9*D(1))) OR ((0.9*SF) >= D(1)) THEN
399: FOR I := 2 UNTIL N DO D(I) := 0;
400: FOR K := 2 UNTIL N DO
401: BEGIN FOR I := 1 UNTIL N DO Y(I) := X(I); SF := FX;
402: ILLC := ILLC OR (KT > 0);
403: L1: KL := K; DF := 0; IF ILLC THEN
404: BEGIN COMMENT: RANDOM STEP TO GET OFF RESOLUTION VALLEY;
405: FOR I := 1 UNTIL N DO
406: BEGIN S := Z(I) := (0.1*LDT + T2*10**KT)*(RANDOM(I)-0.5L);
407: COMMENT: PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM
408: NUMBER UNIFORMLY DISTRIBUTED IN (0, 1) AND
409: THAT ANY INITIALIZATION OF THE RANDOM NUMBER
410: GENERATOR HAS ALREADY BEEN DONE;
411: FOR J := 1 UNTIL N DO X(J) := X(J) + S*V(J,I)
412: END;
413: FX := F(X,N); NF := NF + 1
414: END;
415: FOR K2 := K UNTIL N DO
416: BEGIN SL := FX; S := 0;
417: COMMENT: MINIMIZE ALONG "NON-CONJUGATE" DIRECTIONS;
418: MIN (K2, 2, D(K2), S, FX, FALSE);
419: S := IF ILLC THEN D(K2)*(S + Z(K2))**2 ELSE SL - FX;
420: IF DF < S THEN
421: BEGIN DF := S; KL := K2
422: END
423: END;
424: IF not ILLC AND (DF < ABS( 100*MACHEPS*FX)) THEN
425: BEGIN COMMENT: NO SUCCESS ILLC = FALSE SO TRY ONCE
426: WITH ILLC = TRUE;
427: ILLC := TRUE; GO TO L1
428: END;
429: IF (K = 2) AND (PRIN > 1) THEN VECPRINT ("NEW D", D, N);
430: FOR K2 := 1 UNTIL K - 1 DO
431: BEGIN COMMENT: MINIMIZE ALONG "CONJUGATE" DIRECTIONS;
432: S := 0; MIN (K2, 2, D(K2), S, FX, FALSE)
433: END;
434: F1 := FX; FX := SF; LDS := 0;
435: FOR I := 1 UNTIL N DO
436: BEGIN SL := X(I); X(I) := Y(I); SL := Y(I) := SL - Y(I);
437: LDS := LDS + SL*SL
438: END;
439: LDS := LONGSQRT(LDS); IF LDS > SMALL THEN
440: BEGIN COMMENT: THROW AWAY DIRECTION KL AND MINIMIZE
441: ALONG THE NEW "CONJUGATE" DIRECTION;
442: FOR I := KL - 1 STEP -1 UNTIL K DO
443: BEGIN FOR J := 1 UNTIL N DO V(J,I + 1) := V(J,I);
444: D(I + 1) := D(I)
445: END;
446: D(K) := 0; FOR I := 1 UNTIL N DO V(I,K) := Y(I)/LDS;
447: MIN (K, 4, D(K), LDS, F1, TRUE);
448: IF LDS <= 0 THEN
449: BEGIN LDS := -LDS;
450: FOR I := 1 UNTIL N DO V(I,K) := -V(I,K)
451: END
452: END;
453: LDT := LDFAC*LDT; IF LDT < LDS THEN LDT := LDS;
454: PRINT;
455: T2 := 0; FOR I := 1 UNTIL N DO T2 := T2 + X(I)**2;
456: T2 := M2*LONGSQRT(T2) + T;
457: COMMENT: SEE IF STEP LENGTH EXCEEDS HALF THE TOLERANCE;
458: KT := IF LDT > (0.5*T2) THEN 0 ELSE KT + 1;
459: IF KT > KTM THEN GO TO L2
460: END;
461: COMMENT: TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE STUCK
462: IN A CURVED VALLEY;
463: QUAD;
464: DN := 0; FOR I := 1 UNTIL N DO
465: BEGIN D(I) := 1/LONGSQRT(D(I));
466: IF DN < D(I) THEN DN : = D(I)
467: END;
468: IF PRIN > 3 THEN MATPRINT ("NEW DIRECTIONS", V, N, N);
469: FOR J := 1 UNTIL N DO
470: BEGIN S := D(J)/DN;
471: FOR I := 1 UNTIL N DO V(I,J) := S*V(I,J)
472: END;
473: IF SCBD > 1 THEN
474: BEGIN COMMENT: SCALE AXES TO TRY TO REDUCE CONDITION
475: NUMBER;
476: S := VLARGE; FOR I := 1 UNTIL N DO
477: BEGIN SL := 0; FOR J := 1 UNTIL N DO SL := SL+V(I,J)**2;
478: Z(I) := LONGSQRT(SL);
479: IF Z(I) < M4 THEN Z(I) := M4; IF S > Z(I) THEN S := Z(I)
480: END;
481: FOR I := 1 UNTIL N DO
482: BEGIN SL := S/Z(I); Z(I) := 1/SL; IF Z(I) > SCBD THEN
483: BEGIN SL := 1/SCBD; Z(I) := SCBD
484: END;
485: FOR J := 1 UNTIL N DO V(I,J) := SL*V(I,J)
486: END
487: END;
488: COMMENT: TRANSPOSE V FOR MINFIT LINE BEFORE WAS OMMITTED IN PUBLICATION;
489: FOR I := 2 UNTIL N DO FOR J := 1 UNTIL I - 1 DO
490: BEGIN S := V(I,J); V(I,J) := V(J,I); V(J,I) := S END;
491: COMMENT: FIND THE SINGULAR VALUE DECOMPOSITION OF V, THIS
492: GIVES THE EIGENVALUES AND PRINCIPAL AXES OF THE
493: APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE
494: CONDITION NUMBER;
495: MINFIT (N, MACHEPS, VSMALL, V, D);
496: IF SCBD > 1 THEN
497: BEGIN COMMENT: UNSCALlNG; FOR I := 1 UNTIL N DO
498: BEGIN S := Z(I) ;
499: FOR J := 1 UNTIL N DO V(I,J) := S*V(I,J)
500: END;
501: FOR I := 1 UNTIL N DO
502: BEGIN S := 0; FOR J := 1 UNTIL N DO S := S + V(J,I)**2;
503: S := LONGSQRT(S); D(I) := S*D(I); S := 1/S;
504: FOR J := 1 UNTIL N DO V(J,I) := S*V(J,I)
505: END
506: END;
507: FOR I := 1 UNTIL N DO
508: BEGIN D(I) := IF (DN*D(I)) > LARGE THEN VSMALL ELSE
509: IF (DN*D(I)) < SMALL THEN VLARGE ELSE (DN*D(I))**(-2)
510: END;
511: COMMENT: SORT NEW EIGENVALUES AND EIGENVECTORS;
512: SORT;
513: DMIN := D(N) ; IF DMIN < SMALL THEN DMIN := SMALL;
514: ILLC := (M2*D(1)) > DMIN;
515: IF (PRIN > 1) AND (SCBD > 1) THEN
516: VECPRINT ("SCALE FACTORS", Z, N);
517: IF PRIN > 1 THEN VECPRINT ("EIGENVALUES OF A", D, N);
518: IF PRIN > 3 THEN MATPRINT ("EIGENVECTORS OF A", V, N, N);
519: COMMENT: GO BACK TO MAIN LOOP;
520: GO TO L0;
521: L2: IF PRIN > 0 THEN VECPRINT ("X IS", X, N);
522: FX
523: END ALGOLPRAXIS.
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/small/back.gif){kind=icon}
Return to praxis.alw CVS log ![[TXT]](/icons/small/text.gif){kind=icon}
![[DIR]](/icons/small/dir.gif){kind=icon}
Up to [Local Repository] / imach / src