Annotation of imach/src/praxis.alw, revision 1.1
1.1 ! brouard 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