MODULE Praxis;

IMPORT Wr, Thread, Fmt, Math, Random;
FROM Stdio IMPORT stderr;
  
(* Heuristic numbers:
       If the axes may be badly scaled (wich is to be avoided if possible).
       Then set scbd = 10, otherwise set scbd = 1. If the problem is known
       to be ill-conditioned, set illc = TRUE, otherwise set illc = FALSE.
       ktm is the number of iterations without improviment before the 
       algorithm terminates. ktm = 4 is very cautious: usually ktm = 1 is
       satisfactory.
*)

PROCEDURE Minimize(
  f: GoalFunc;
  VAR x: ARRAY OF REAL;   (* In: starting point, Out: best point found. *)
  VAR fx: REAL;           (* Out: value of "f" at optimum. *)
  t0: REAL := 1.0e-6;     (* tolerance: norm(x-x0) < t0 + sqrt(machep)*norm(x) *)
  scbd: REAL := 1.0;      (* scbd = 10.0 if the axes may be badly scaled, is to be avoided if possible.*)
  illc: BOOLEAN := FALSE; (* TRUE if the problem is known to be ill-conditioned. *)
  ktm: CARDINAL := 1;     (* number of iterations without improvement before the algorithms terminates *)
  prin: CARDINAL := 0;    (* prin = 0, nothing is printed. prin = [0, 1, 2, 3, 4] *)
  maxCalls: CARDINAL;     (* Maximum number of calls to "f". *)
  machep: REAL := 1.0e-15; (* Minimum number size. *)
  maxStep: REAL := 1.0e+30; (* Maximum step size *)
  report: ReportProc := NIL;       (* Called when a new minimum is found. *)
  normalize: NormalizeProc := NIL; (* Called to modify point before applying "f". *)
): REAL  =  

  VAR Div0: REAL := machep; (* Avoid division for zero *)
  BEGIN
    WITH
      N = NUMBER(x),
      n = N - 1,
      v = NEW(REF ARRAY OF ARRAY OF REAL, N,N)^,
      q0 = NEW(REF ARRAY OF REAL, N)^,
      q1 = NEW(REF ARRAY OF REAL, N)^,
      d = NEW(REF ARRAY OF REAL, N)^,
      y = NEW(REF ARRAY OF REAL, N)^,
      z = NEW(REF ARRAY OF REAL, N)^,
      xto = NEW(REF ARRAY OF REAL, N)^ (* save the best solution at moment *)
    DO
 
      VAR qd0, qd1, qf1: REAL := 0.0; (* Used in procedure Min and/or Quad *)
          dmin, ldt: REAL := 0.0; (* Used in procedure Min, too *)
          small, vsmall, large, vlarge, s, sl, sf, dn, ft, f1, t, h,
          df, m2, m4, ldfac, lds, dni, t2: REAL := 0.0;
          calls, i, im1, kl, klmk, km1, nf, nl: CARDINAL := 0;
          kt: INTEGER := 0;
          boole: BOOLEAN := FALSE;
      
      PROCEDURE CallF (VAR z: ARRAY OF REAL; VAR u: REAL) =
        (* Calls normalize on "z", then evaluates "f(z, u, calls)". *)
        BEGIN
          IF normalize # NIL THEN normalize(z) END;
          INC(calls);
          f(z, u, calls);
          IF u < ft THEN 
            ft := u;           
            FOR i := 0 TO N-1  DO xto[i] := z[i] END
          END;
        END CallF;
   
      PROCEDURE MinFit() =
    (* An improved version of MinFit (see Golub and Reinsch, 1969) restricted*)
    (* to m = n, p=o. 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. In this implementation   *)
    (* ab = v and q = d.                                                     *)

    (* Householder's reduction to bidiagonal form.                           *)
    
      VAR boole, boole2: BOOLEAN := FALSE;
          d := NEW(REF ARRAY OF REAL, N);
          e := NEW(REF ARRAY OF REAL, N); 
          c, eps, f2, g, h, s, temp, x2, y2, z: REAL := 0.0;
          i, k, kt, l, l2, lp1: CARDINAL := 0;
          ll2: INTEGER :=  0;
      <* FATAL Wr.Failure, Thread.Alerted *>       
      BEGIN
        IF (n # 0) THEN
          eps := machep;
          g := 0.0;
          x2 := 0.0;
          FOR i := 0 TO N-1 DO
            e[i] := g;
            s := 0.0;
            l := i + 1;
            FOR j := 0 TO N-1 DO s := s + v[j][i]*v[j][i] END;
            g := 0.0;
            IF (s >= vsmall) THEN 
              f2 := v[i][i];
              g := FLOAT(Math.sqrt(FLOAT(s, LONGREAL))); 
              IF (f2 >= 0.0) THEN g := -g END;
              h := f2*g - s;
              v[i][i] := f2 - g;
              IF (l <= n) THEN
                FOR j := l TO N-1 DO
                  f2 := 0.0;
                  FOR k := i TO N-1 DO f2 := f2 + v[k][i]*v[k][j] END;
                  f2 := f2/(h+Div0);
                  FOR k := i TO N-1 DO v[k][j] := v[k][j] + f2*v[k][i] END;
                END;
              END;
            END;
            d[i] := g;
            s := 0.0;
            IF (i # n) THEN
              FOR j := l TO N-1 DO s := s + v[i][j]*v[i][j] END;
            END;
            g := 0.0;
            IF (s >= vsmall) THEN 
              IF (i # n) THEN f2 := v[i][i+1]  END;
              g := FLOAT(Math.sqrt(FLOAT(s, LONGREAL)));
              IF (f2 >= 0.0) THEN g := -g END;
              h := f2*g -s;
              IF (i # n) THEN 
                v[i][i+1] := f2 - g;
                FOR j := l TO N-1 DO e[j] := v[i][j]/(h+Div0) END;
                FOR j := l TO N-1 DO
                  s := 0.0;
                  FOR k := l TO N-1 DO s := s + v[j][k]*v[i][k] END;
                  FOR k := l TO N-1 DO v[j][k] := v[j][k] + s*e[k] END;
                END;
              END;
            END;
            y2 := ABS(d[i]) + ABS(e[i]);
            IF (y2 > x2) THEN x2 := y2 END;
          END;
    (* accumulation of right-hand transformations   *)
          v[n][n] := 1.0;
          g := e[n];
          l := n;
          FOR ii := 1 TO N-1 DO
            i := n - ii ;
            IF (g # 0.0) THEN 
              h := v[i][i+1]*g;
              FOR j := l TO N-1 DO v[j][i] := v[i][j]/(h+Div0) END;
              FOR j := l TO N-1 DO 
                s := 0.0;
                FOR k := l TO N-1 DO s := s+v[i][k]*v[k][j] END;
                FOR k := l TO N-1 DO v[k][j] := v[k][j] + s*v[k][i] END;
              END;
            END;
            FOR j := l TO N-1 DO
              v[i][j] := 0.0;
              v[j][i] := 0.0;
            END;
            v[i][i] := 1.0;
            g := e[i];
            l := i;
          END;
    (* Diagonalization of the bidiagonal form *)        
          eps := eps*x2;
          FOR kk := 0 TO N-1 DO
            k := n - kk;
            kt := 0;
            boole := FALSE;
            boole2 := FALSE;
            REPEAT
              kt := kt + 1;
              IF (kt > 30) THEN 
                e[k] := 0.0;
                Wr.PutText(stderr, " qr failed");
                Wr.PutText(stderr, "\n");
              END;
              ll2 := -1;
              REPEAT
                ll2 := ll2 + 1;
                l2 := k - ll2;
                l := l2;
                boole := (ABS(e[l]) <= eps);
                IF (l # 0) THEN
                  boole2 := (ABS(d[l-1]) <= eps);
                END;
              UNTIL ( (ll2 = k) OR ( boole) OR (boole2) );
              boole2 := FALSE;

              IF NOT(boole) THEN 
    (* Cancellation of e[l] if l>0 *)
                c := 0.0;
                s := 1.0;
                FOR i := l TO k DO
                  f2 := s*e[i];
                  e[i] := c*e[i];
                  IF (ABS(f2) > eps) THEN 
                    g := d[i];
    (* q[i] = h = sqrt(g*g + f2*f2)  *)
                    IF (ABS(f2) >= ABS(g)) THEN 
                      IF (ABS(f2) > 0.0) THEN
                        h := ABS(f2)*( FLOAT(Math.sqrt(FLOAT( 
                                          (1.0 + (g*g)/(f2*f2+Div0)), LONGREAL))) );
                      ELSE
                        h := 0.0;
                      END;
                    ELSE
                      h := ABS(g)*( FLOAT(Math.sqrt(FLOAT(
                                          (1.0 + (f2*f2)/(g*g+Div0)), LONGREAL))) );
                    END;
                    d[i] := h;
                    IF h = 0.0 THEN
                      g := 1.0;
                      h := 1.0;
                    END;
                    c := g/(h+Div0);
                    s := -f2/(h+Div0);
                  END;
                END;
              END;
              boole := FALSE;            
    (* Test for convergence *)
              z := d[k];
              boole := (l # k);
              IF (boole) THEN
    (* Shift from bottom 2*2 minor *)
                x2 := d[l];
                y2 := d[k-1];
                g := e[k-1];
                h := e[k];
                f2 := ((y2 - z)*(y2 + z) + (g - h)*(g + h))/(2.0*h*y2+Div0);
                g := FLOAT(Math.sqrt(FLOAT( (1.0 + (f2*f2)), LONGREAL)));
                temp := f2 - g;
                IF (f2 >= 0.0)  THEN temp := f2 + g END;
                f2 := ((x2 - z)*(x2 + z) + h*(y2/(temp+Div0) - h))/(x2+Div0);
    (* Next qr transformation *)
                c := 1.0;
                s := 1.0;
                lp1 := l + 1;
                IF (lp1 <= k) THEN
                  FOR i := lp1 TO k DO
                    g := e[i];
                    y2 := d[i];
                    h := s*g;
                    g := g*c;
                    IF (ABS(f2) >= ABS(h)) THEN 
                      IF (ABS(f2) > 0.0) THEN
                        z := ABS(f2)*( FLOAT(Math.sqrt(FLOAT( 
                                          (1.0 + (h*h)/(f2*f2+Div0)), LONGREAL))) );
                      ELSE
                        z := 0.0;
                      END;
                    ELSE
                      z := ABS(h)*( FLOAT(Math.sqrt(FLOAT( 
                                          (1.0 + (f2*f2)/(h*h+Div0)), LONGREAL))) );
                    END;
                    e[i-1] := z;
                    IF (z = 0.0) THEN
                      f2 := 1.0;
                      z := 1.0;
                    END;
                    c := f2/(z+Div0);
                    s := h/(z+Div0);
                    f2 := x2*c + g*s;
                    g := -x2*s + g*c;
                    h := y2*s;
                    y2 := y2*c;
                    FOR j := 0 TO N-1 DO
                      x2 := v[j][i-1];
                      z := v[j][i];
                      v[j][i-1] := x2*c + z*s;
                      v[j][i] := -x2*s + z*c
                    END;
                    IF (ABS(f2) >= ABS(h)) THEN 
                      IF (ABS(f2) > 0.0) THEN
                        z := ABS(f2)*( FLOAT(Math.sqrt(FLOAT( 
                                          (1.0 + (h*h)/(f2*f2+Div0)), LONGREAL))) );
                      ELSE
                        z := 0.0;
                      END;
                    ELSE
                      z := ABS(h)*( FLOAT(Math.sqrt(FLOAT( 
                                          (1.0 + (f2*f2)/(h*h+Div0)), LONGREAL))) );
                    END;
                    d[i-1] := z;
                    IF (z = 0.0) THEN
                      f2 := 1.0;
                      z := 1.0;
                    END;
                    c := f2/(z+Div0);
                    s := h/(z+Div0);
                    f2 := c*g + s*y2;
                    x2 := -s*g + c*y2;
                  END;
                END;
                e[l] := 0.0;
                e[k] := f2;
                d[k] := x2;
              END;
            UNTIL NOT(boole);
    (* Convergence: d[k] is made non-negative *)
            IF (z < 0.0) THEN 
              d[k] := -z;
              FOR j := 0 TO N-1 DO v[j][k] := -v[j][k] END;
            END;
          END;
        ELSE
          d[0] := v[0][0];
          v[0][0] := 0.0;
        END;
        Print(x, "*** MinFit ***");        
      END MinFit;

      PROCEDURE Min(j: INTEGER; nits: CARDINAL; VAR d2, 
                    x1: REAL; f1: REAL; fk: BOOLEAN)  =
    (*   The subroutine min minimizes f from x in the direction v[*][j]      *)
    (*   unless j is less than 1, when a quadratic search is made in the     *)
    (*   plane defined by q0,q1,x.                                           *)
    (*   d2 is either zero or an approximation to half f".                   *)
    (*   On entry, x1 is an estimate of the distance from x to the minimum   *)
    (*   along v[*][j] (or, if j=0, a curve).  On return, x1 is the distance *)
    (*   found.                                                              *)
    (*   If fk=.true., then f1 is flin(x1).  otherwise x1 and f1 are ignored *)
    (*   on entry unless final fx is greater than f1.                        *)
    (*   nits controls the number of times an attempt will be made to halve  *)
    (*   the interval.                                                       *)

        PROCEDURE Flin(j: INTEGER; VAR l: REAL): REAL  =
        VAR xt := NEW(REF ARRAY OF REAL, N);
            qa, qb, qc: REAL := 0.0;
        BEGIN
          IF (j # -1) THEN 
    (* The search is linear *)      
            FOR i := 0 TO N-1 DO 
              xt[i] := x[i] + l*v[i][j];
            END;
          ELSE
    (* The search is along a parabolic space curve *)
            qa:= (l*(l - qd1))/(qd0*(qd0 + qd1)+Div0);
            qb := ((l + qd0)*(qd1 - l))/(qd0*qd1+Div0);
            qc := (l*(l + qd0))/(qd1*(qd0 + qd1)+Div0);
            FOR i := 0 TO N-1 DO
              xt[i] := (qa*q0[i] + qb*x[i]) + qc*q1[i];
            END;
          END;  
          
    (* The function evaluation counter nf is incremented *)
          nf := nf + 1;
          Print(xt^, "Print I ** Flin");
          CallF(xt^, fx);
          Print(xt^, "Print E ** Flin");
          RETURN fx;
        END Flin; 

      VAR dz, boole, boole2: BOOLEAN := FALSE;    
          k: CARDINAL := 0;
          d1, f0, f2, fm, s, sf1, sx1, t2, temp, x2, xm: REAL := 0.0;
      BEGIN
        f1 := f1;
        sf1 := f1;
        sx1 := x1;
        k := 0;
        xm := 0.0;
        fm := fx;
        f0 := fx;
        dz := (d2 < machep);
    (* Find the step size  *)
        s := 0.0;
        FOR i := 0 TO N-1 DO s := s + x[i]*x[i] END;
        s := FLOAT(Math.sqrt(FLOAT(s, LONGREAL)));
        temp := d2;
        IF (dz) THEN temp := dmin END;
        t2 :=  m4*( FLOAT(Math.sqrt(FLOAT(ABS(fx)/(temp+Div0) + s*ldt, LONGREAL))) ) + m2*ldt;
        s := m4*s + t;
        IF ( (dz) AND (t2 > s) ) THEN t2 := s END;
        t2 := MAX(t2, small);
        t2 := MIN(t2, 0.01*h);
        IF NOT(NOT(fk) OR (f1 > fm)) THEN
          xm := x1;
          fm := f1;
        END;
        IF NOT( (fk) AND (ABS(x1) >= t2) )  THEN
          temp := 1.0;
          IF (x1 < 0.0) THEN temp := -1.0 END;
          x1 := temp*t2;
          f1 := Flin(j, x1);
        END;
        IF (f1 <= fm) THEN 
          xm := x1;
          fm := f1;
        END;
        REPEAT
          IF (dz) THEN 
    (* Evaluate Flin at another point and estimate the second derivative *)
            x2 := -x1;
            IF (f0 >= f1) THEN x2 := 2.0*x1 END;
            f2 := Flin(j, x2);
            IF (f2 <= fm) THEN 
              xm := x2;
              fm := f2;
            END;
            d2:= (x2*(f1 - f0)-x1*(f2 - f0))/((x1*x2)*(x1 - x2)+Div0); 
          END;        
    (* Estimate the first derivative at 0 *)
          d1 := (f1 - f0)/(x1+Div0) - x1*d2;
          dz := TRUE;
    (* Predict the minimum *)
          IF (d2 <= small) THEN 
            x2 := h;
            IF (d1 >= 0.0) THEN x2 := -x2 END;
          ELSE
            x2 := (-0.5*d1)/(d2+Div0);
          END;
          IF (ABS(x2) > h) THEN
            IF x2 <= 0.0 THEN 
              x2 := -h;
            ELSE
              x2 := h;
            END;
          END;
    (* Evaluate f at the predicted minimum *)
          boole := FALSE;
          boole2 := FALSE;
          WHILE NOT(boole2) DO 
            boole2 := TRUE;
            f2 := Flin(j, x2);  
            boole := ( (k >= nits) OR (f2 <= f0) );
    (* No success. So try again *)
            IF NOT(boole) THEN 
              k := k + 1;
              boole2 := ( (f0 < f1) AND ((x1*x2) > 0.0) );
              IF NOT(boole2) THEN 
                x2 := 0.5*x2;
              END;
            END;
          END;
        UNTIL (boole);
    (* Increment the one-dimensional search counter *)      
        nl := nl + 1;
        IF (f2 > fm)  THEN
          x2 := xm;
        ELSE
          fm := f2;
        END;
    (* Get a new estimate fo the second derivative *) 
        IF ( ABS(x2*(x2 - x1)) > small) THEN 
          d2:= (x2*(f1-f0) - x1*(fm-f0))/((x1*x2)*(x1-x2)+Div0);  
        ELSE
          IF (k > 0) THEN d2 := 0.0 END;
        END;
        IF (d2 <= small) THEN d2 := small END;
        x1 := x2;
        fx := fm;
        IF (sf1 < fx) THEN 
          fx := sf1;
          x1 := sx1;
        END;
    (* Update x for linear but not parabolic search *)
        IF (j # -1) THEN 
          FOR i := 0 TO N-1 DO x[i] := x[i] + x1*v[i][j] END;
        END;
        Print(x, "Print  ** Min");
      END Min;

      PROCEDURE Sort() =
    (* Sorts the elements of d[n] into descending order and moves the *)
    (* corresponding columns od v[n][n]. m is the row dimension of v  *)
    (* as declared in the calling program.                            *)
    
      VAR k, ip1, nm1: CARDINAL :=0;
          s: REAL := 0.0;
      BEGIN  
        IF (n # 0) THEN 
          nm1 := n - 1;
          FOR i := 0 TO nm1 DO
            k := i;
            s := d[i];
            ip1 := i + 1;
            FOR j := ip1 TO N-1 DO
              IF (d[j] > s) THEN 
                k := j;
                s := d[j]
              END;
            END;
            IF (k > i) THEN
              d[k] := d[i];
              d[i] := s;
              FOR j := 0 TO N-1 DO
                s := v[j][i];
                v[j][i] := v[j][k];
                v[j][k] := s
              END;
            END;
          END;
        END;
      END Sort;

      PROCEDURE Quad() =
      VAR l, qa, qb, qc, qf1, s: REAL := 0.0;
    (* Quad looks for the minimum of f along a curve defined by q0, q1, x *)
      BEGIN
        s := fx;
        fx := qf1;
        qf1 := s;
        FOR i := 0 TO N-1 DO
          s := x[i];
          l := q1[i];
          x[i] := l;
          q1[i] := s;
          qd1 := qd1 + (s-l)*(s-l);
        END;
        qd1 := FLOAT(Math.sqrt(FLOAT(qd1, LONGREAL)));
        l := qd1;
        s := 0.0;
        IF NOT((qd0 <= 0.0) OR (qd1 <= 0.0) OR (nl < 3*n*n)) THEN
          Min(-1, 2, s, l, qf1, TRUE);
          
          IF calls > maxCalls THEN  
            (* The number of calls of func exceded maxCalls *)
            RETURN
          END;
          
          qa := (l*(l-qd1))/(qd0*(qd0+qd1)+Div0);
          qb := ((l+qd0)*(qd1-l))/(qd0*qd1+Div0);
          qc := (l*(l+qd0))/(qd1*(qd0+qd1)+Div0);
        ELSE
          fx := qf1;
          qa := 0.0;
          qb := qa;
          qc := 1.0;
        END;
        qd0 := qd1; 
        FOR i := 0 TO N-1 DO
          s := q0[i];
          q0[i] := x[i];
          x[i] := (qa*s + qb*x[i]) + qc*q1[i];
        END;
      END Quad;

      PROCEDURE VCPrnt(option: CARDINAL; e: ARRAY OF REAL) =
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN
        CASE option OF 
        |1 =>
             Wr.PutText(stderr, "The second diference array d[*] is:" & "\n");
             FOR i := 0 TO N-1 BY 3 DO
               Wr.PutText(stderr, Fmt.Real(e[i]) &" "& Fmt.Real(e[i+1]) &" "
                          &  Fmt.Real(e[i+2]) &"\n");
             END;
        |2 =>
             Wr.PutText(stderr, "The scale factors are:" & "\n");
             FOR i := 0 TO N-1 BY 3 DO
               Wr.PutText(stderr, Fmt.Real(e[i]) &" "& Fmt.Real(e[i+1]) &" "
                          &  Fmt.Real(e[i+2]) &"\n");
             END;
        |3 =>
             Wr.PutText(stderr, "The approximating quadratic form has the principal values:"  & "\n");
             FOR i := 0 TO N-1 BY 3 DO
               Wr.PutText(stderr, Fmt.Real(e[i]) &" "& Fmt.Real(e[i+1]) &" "
                          &  Fmt.Real(e[i+2]) &"\n");
             END;
        |4 =>
             Wr.PutText(stderr, "x is:" & "\n");
             FOR i := 0 TO N-1 BY 3 DO
               Wr.PutText(stderr, Fmt.Real(e[i]) &" "& Fmt.Real(e[i+1]) &" "
                          &  Fmt.Real(e[i+2]) &"\n");
             END;
        END; 
      END VCPrnt;

      PROCEDURE Print(X: ARRAY OF REAL; texto: TEXT) =
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN
         IF NOT( (n >= 2) AND ( prin <= 2) ) THEN
           Wr.PutText(stderr, texto & "\n");
           Wr.PutText(stderr, "After " & Fmt.Int(nl) & " linear searches, the function has been evaluated " & Fmt.Int(nf) & " times. The smallest value found is f[x] =  " & Fmt.Real(fx) & "\n");
           
           FOR i := 0 TO N-1 BY 3 DO
             Wr.PutText(stderr, Fmt.Real(X[i]) &" "& Fmt.Real(X[i+1]) &" "
                        &  Fmt.Real(X[i+2]) &"\n");
           END;
         END;
      END Print;

      PROCEDURE MAPrnt(option: CARDINAL) =
    (* The subroutine maprnt prints the columns of the nxn matrix v    *)
    (* with a heading as specified by option.                          *)
    (* m is the row dimension of v as declared in the calling program  *)

      VAR boole: BOOLEAN := FALSE;
          low, upp: CARDINAL := 0;
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN
        low := 1;
        upp := 5;
        CASE option OF 
        |1 =>
             Wr.PutText(stderr, " The NEW directions are:" & "\n");
        |2 =>
             Wr.PutText(stderr, " AND the principal axes:" & "\n");
        END;
        REPEAT
          IF (n < upp) THEN upp := n END;
          FOR i := 0 TO N-1 DO
            FOR j := low TO upp BY 3 DO
              Wr.PutText(stderr, Fmt.Real(v[i][j]) &" "& Fmt.Real(v[i][j+1])
                         &" " &  Fmt.Real(v[i][j+2]) &"\n");
            END;
          END;
          low := low + 5;
          boole := (n >= low);
          IF (boole) THEN 
            upp := upp + 5;
          END;
        UNTIL NOT(boole);
      END MAPrnt;
    
      BEGIN 
  (* Initialization - Machine Dependent Numbers *)
        small := machep*machep;
        vsmall := small*small;
        large := 1.0/small;
        vlarge := 1.0/vsmall;
        m2 := FLOAT( Math.sqrt(FLOAT(machep, LONGREAL)) );
        m4 := FLOAT( Math.sqrt(FLOAT(m2, LONGREAL)) );

        ldfac := 0.01;
        IF (illc) THEN ldfac := 0.1 END;
        kt := 0;
        nl := 0;
        nf := 1;
        CallF(x, fx);
        qf1 := fx;
        ft := fx;
        t := small + ABS(t0);
        t2 := t;
        dmin := small;
        h := maxStep;
        IF (h < 100.0*t0) THEN h := 100.0*t END;
        ldt := h;

  (* The first set of search directions "v" is the identity matrix *)
        FOR i := 0 TO N-1  DO 
          FOR j := 0 TO N-1 DO
            v[i][j] := 0.0;
          END;
          v[i][i] := 1.0;
          d[i] := 0.0;
        END;
        qd0 := 0.0;
        FOR i := 0 TO N-1 DO
          q0[i] := x[i];
          q1[i] := x[i];
        END;

        IF (prin > 0) THEN  Print(x, "Main - Inicio") END;
  (* ***************** The main loop starts here *****************  *)
        WHILE calls < maxCalls DO
          sf := d[0];
          d[0] := 0.0;
          s := 0.0;
  (* Minimize along the FIRST direction v[*][0] *)          
          Min(0, 2, d[0], s, fx, FALSE);
  
          IF calls > maxCalls THEN  
            (* The number of calls of func exceded maxCalls *)
            RETURN ft
          END;
          
          IF (s <= 0.0) THEN 
            FOR i := 0 TO N-1 DO v[i][0] := -v[i][0] END;
          END;
          IF NOT( (sf > 0.9*d[0]) AND (0.9*sf < d[0]) ) THEN
            FOR i := 1 TO N-1 DO d[i] := 0.0 END;
          END;
  (* ...... The inner loop starts here *)        
          FOR k := 1 TO N-1 DO (* *** Begin big FOR *** *) 
            Print(x, "Inner Loop Starts");
            FOR i := 0 TO N-1 DO 
              y[i] := x[i];
              xto[i] := x[i]
            END;
            sf := fx;
            IF (kt > 0) THEN illc := TRUE END;
            REPEAT
              kl := k;
              df := 0.0;
  (* A random step follows (to avoid resolution valleys), praxis assumes that *)
  (* random returns a random number uniformly distributed in (0,1).           *)
              IF illc THEN 
                FOR i := 0 TO N-1 DO 
                  s := 0.1*ldt + t2 * FLOAT(
                                       Math.pow(10.0D0, FLOAT(kt, LONGREAL)) ) 
                                     * (Random.Real() - 0.5);
                  z[i] := s;
                  FOR j := 0 TO N-1 DO
                    x[j] := x[j] + s*v[j][i];
                  END;
                END;
                Print(x, "Random I");
                CallF(x, fx);
                Print(x, "Random E");
                nf := nf + 1;
              END;              
  (* Minimize along the "non-conjugate" directions v[*][k], ..., v[*][n] *)
              FOR k2 := k TO N-1 DO
                sl := fx;
                s := 0.0;               
                Min(k2, 2, d[k2], s, fx, FALSE);
                
                IF calls > maxCalls THEN  
                (* The number of calls of func exceded maxCalls *)
                  RETURN ft
                END;
                 
                IF illc THEN
                  s := d[k2]*( (s+z[k2])*(s+z[k2]) );
                ELSE
                  s := sl - fx;
                END;
                IF (df <= s) THEN
                  df := s;
                  kl := k2;
                END;                  
              END; 
              boole := ( illc OR ( df >= ABS(100.0*machep*fx) ) );
              IF NOT(boole) THEN illc := TRUE END;
            UNTIL (boole);
            IF (k = 2) AND (prin > 1) THEN VCPrnt(1, d) END;
  (* Minimize along the "conjugate" directions v[*][0], ..., v[*][k-1] *)
            km1 := k - 1;            
            FOR k2 := 0 TO km1 DO 
              s := 0.0;
              Min(k2, 2, d[k2], s, fx, FALSE);

              IF calls > maxCalls THEN  
              (* The number of calls of func exceded maxCalls *)
                RETURN ft
              END;
            END;
            f1 := fx;
            fx := sf;
            lds := 0.0;
            FOR i := 0 TO N-1 DO
              sl := x[i];
              x[i] := y[i];
              sl := sl - y[i];
              y[i] := sl;
              lds := lds + sl*sl;
            END;
            lds := FLOAT( Math.sqrt(FLOAT(lds, LONGREAL)) );
            IF (lds > small) THEN
  (* Discard direction v[*][kl]. If no random step was taken, v[*][kl] is the *)
  (* "non-conjugate" direction alonh which the greatest improvement was made. *)
              klmk := kl - k;
              IF (klmk >= 0) THEN
                FOR ii := 0 TO klmk DO
                  i := kl - ii - 1;
                  FOR j := 0 TO N-1 DO v[j][i+1] := v[j][i] END;
                  d[i+1] := d[i];
                END;
              END;
              d[k] := 0.0;
              FOR i := 0 TO N-1 DO v[i][k] := y[i]/(lds+Div0) END;
  (* Minimize along the new "conjugate" direction v[*][k], which is the  *)
  (* normalized vector: (new x) - (old x)                                *)
              Min(k, 4, d[k], lds, f1, TRUE);
              
              IF calls > maxCalls THEN  
              (* The number of calls of func exceded maxCalls *)
                RETURN ft
              END;
              
              IF (lds <= 0.0) THEN
                lds := -lds;
                FOR i := 0 TO N-1 DO v[i][k] := -v[i][k] END;
              END;
            END;
            ldt := ldfac*ldt;
            IF (ldt < lds) THEN ldt := lds END;
            IF (prin > 0) THEN Print(x, "Min. along new conjg. dir. v[*][k]") END;
            t2 := 0.0;
            FOR i := 0 TO N-1 DO t2 := t2 + x[i]*x[i] END;
            t2 := m2*FLOAT( Math.sqrt(FLOAT(t2, LONGREAL)) ) + t;
  (* See whether the length of step taken since starting the inner loop *)
  (* exceeds half the tolerance.                                        *)
            IF (ldt > (0.5*t2)) THEN kt := -1 END;
            kt := kt + 1;
            IF (kt > ktm) THEN
              IF (prin > 0) THEN VCPrnt(4, x) END;
              RETURN fx;
            END;
          END;      (* *** End big FOR **** *)
  (* ...... The inner loop ends here  *)
  (*        Try quadratic extrapolation in case we are in a curved valley *)
          Quad();
          
          IF calls > maxCalls THEN  
          (* The number of calls of func exceded maxCalls *)
            RETURN ft
          END;
          
          dn := 0.0;
          FOR i := 0 TO N-1 DO
            d[i] := 1.0/( FLOAT(Math.sqrt(FLOAT(d[i], LONGREAL))) + Div0 );
            IF (dn < d[i]) THEN dn := d[i] END;
          END;
          IF (prin > 3) THEN MAPrnt(1) END;
          FOR j := 0 TO N-1 DO
            s := d[j]/(dn+Div0);
            FOR i:= 0 TO N-1 DO v[i][j] := s*v[i][j] END;
          END;
  (* Scale the axes to try to reduce the condition NUMBER  *)
          IF (scbd > 1.0) THEN 
            s := vlarge;
            FOR i := 0 TO N-1 DO
              sl := 0.0;
              FOR j := 0 TO N-1 DO sl := sl + v[i][j]*v[i][j] END;
              z[i] := FLOAT(Math.sqrt(FLOAT(sl, LONGREAL)));
              IF (z[i] < m4) THEN z[i] := m4 END;
              IF (s > z[i]) THEN s := z[i] END;
            END;
            FOR i := 0 TO N-1 DO
              sl := s/(z[i]+Div0);
              z[i] := 1.0/(sl+Div0);
              IF (z[i] > scbd) THEN 
                sl := 1.0/(scbd+Div0);
                z[i] := scbd;
              END;
              FOR j := 0 TO N-1 DO v[i][j] := sl*v[i][j] END;
            END;
          END;
  (* Calculate a new set of orthogonal directions before repeating the main *)
  (* loop. First, transpose v for minfit:                                   *)
          FOR i := 1 TO N-1 DO
            im1 :=  i - 1;
            FOR j := 0 TO im1 DO
              s := v[i][j];
              v[i][j] := v[j][i];
              v[j][i] := s;
            END;
          END;
  (* Call MinFit to find the singular value decomposition of v, this gives the *)
  (* principal values and principal directions of the aproximating quadratic   *)
  (* form without squaring the condition NUMBER                                *)
          MinFit();
  (* Unscale the axes *)
          IF (scbd > 1.0) THEN
            FOR i := 0 TO N-1 DO
              s := z[i];
              FOR j := 0 TO N-1 DO v[i][j] := s*v[i][j] END;
            END;
            FOR i:= 0 TO N-1 DO 
              s := 0.0;
              FOR j := 0 TO N-1 DO s := s + v[j][i]*v[j][i] END;
              s := FLOAT(Math.sqrt(FLOAT(s, LONGREAL)));
              d[i] := s*d[i];
              s := 1.0/(s+Div0);
              FOR j := 0 TO N-1 DO v[j][i] := s*v[j][i] END;
            END;
          END;
          FOR i := 0 TO N-1 DO
            dni := dn*d[i];
            IF (dni <= large) THEN 
              IF (dni >= small) THEN
                d[i] := 1.0/(dni*dni+Div0);
              ELSE
                d[i] := vlarge;
              END;
            ELSE
              d[i] := vsmall;
            END;
          END;
  (* Sort the eigenvalues AND eigenvectors  *)
          Sort();
          dmin := d[n];
          IF (dmin < small) THEN dmin := small END;
          illc := FALSE;
          IF (m2*d[0] > dmin) THEN illc := TRUE END;
          IF (prin > 1) AND (scbd > 1.0) THEN VCPrnt(2, z) END;
          IF (prin > 1) THEN VCPrnt(3, d) END;
          IF (prin > 3) THEN MAPrnt(2) END;
  (* ***************** The main loop ends here *****************  *)      
        END;
      END; 
    END;
    RETURN fx;  
  END Minimize;

BEGIN

END Praxis.