Return to praxis.alw CVS log

Up to [Local Repository] / imach / src

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.

  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.