{Ukazka pomalosti hledani minimaci primym pouzitim gradientu
}
CONST
   np =  10;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
FUNCTION fnc(VAR p: RealArrayNP): real;
BEGIN
     FNC:= p[1]*(101*p[1]-198*p[2]+4)+p[2]*(101*p[2]+4) +5;
END;
PROCEDURE dfnc(VAR p,g: RealArrayNP);
BEGIN
     g[1]:= 202*p[1]-198*p[2] + 4;
     g[2]:= 202*p[2]-198*p[1] +4;
END;

PROCEDURE frprmn(VAR p: RealArrayNP;
                     n: integer;
                  ftol: real;
              VAR iter: integer;
              VAR fret: real);
{Hleda minimum funkce FNC(X) N promennych (X je vektor dimenze N)
 ktera ma gradient dany funkci DFNC pomoci Fletcher-Reeves-Polak-Ribiere
 metody konjugovanych gradientu. Pocatecni bod je P[1..N].
 Relativni presnost minima je FTOL. Na konci P je nalezena hodnota minima,
 ITER je pocet provedenych iteraci, FRET je funkcni hodnota v minimu.
 Procedura pouziva LINMIN pro hledani minima ve smeru.
 ITMAX je maximalni pocet iteraci, EPS je absolutni presnost, pokud
 je minimum funkce velmi blizke 0. }

VAR
   LinminNcom: integer;
   LinminPcom,LinminXicom: RealArrayNP;

PROCEDURE linmin(VAR p,xi: RealArrayNP;
                        n: integer;
                 VAR fret: real);
{V N promennych je dan bod P[1..N] a N dimensionalni smer XI[1..N]
 hleda minimum ve smeru, ve vysledku je P minimum ve smeru XI,
 FRET = FNC(P). Pouziva MNBRAK a BRENT, TOL je hodnopa relativni
 presnosti pro BRENT.}
CONST
   tol = 1.0e-4;
VAR
   j: integer;
   xx,xmin,fx,fb,fa,bx,ax: real;

FUNCTION f1dim(x: real): real;
{Musi doprovazet LINMIN. Prevadi FNC funkci v N promennych na
 jednorozmernou funkci F1DIM.}
VAR
   j: integer;
   xt: ^RealArrayNP;
BEGIN
   new(xt);
   FOR j := 1 TO LinminNcom DO
      xt^[j] := LinminPcom[j]+x*LinminXicom[j];
   f1dim := fnc(xt^);
   dispose(xt)
END;

FUNCTION func(x: real): real;
BEGIN
   func := f1dim(x)
END;


FUNCTION brent(ax,bx,cx,tol: real;
                   VAR xmin: real): real;
{Hleda Brentovou metodou
 minimum zadane funkce FUNC, ktere je jiz dopredu ohraniceno
 hodnotami AX,BX,CX, kde BX lezi mezi AX a CX a funkce je v BX minimalni
 z dane trojice. TOL zadava relativni presnost hledani minima,
 XMIN je na vystupu nalezena hodnota minima, hodnota funkce v minimu
 je BRENT.}
LABEL 99;
CONST
   itmax = 100;
   cgold = 0.3819660;
   zeps = 1.0e-10;
VAR
   a,b,d,e,etemp: real;
   fu,fv,fw,fx: real;
   iter: integer;
   p,q,r,tol1,tol2: real;
   u,v,w,x,xm: real;

FUNCTION sign(a,b: real): real;
BEGIN
   IF b >= 0.0 THEN sign := abs(a) ELSE sign := -abs(a)
END;

BEGIN
   IF ax < cx THEN a := ax ELSE a := cx;
   IF ax > cx THEN b := ax ELSE b := cx;
   v := bx;
   w := v;
   x := v;
   e := 0.0;
   fx := func(x);
   fv := fx;
   fw := fx;
   FOR iter := 1 TO itmax DO BEGIN
      xm := 0.5*(a+b);
      tol1 := tol*abs(x)+zeps;
      tol2 := 2.0*tol1;
      IF abs(x-xm) <= tol2-0.5*(b-a) THEN GOTO 99;
      IF abs(e) > tol1 THEN BEGIN
         r := (x-w)*(fx-fv);
         q := (x-v)*(fx-fw);
         p := (x-v)*q-(x-w)*r;
         q := 2.0*(q-r);
         IF q > 0.0 THEN p := -p;
         q := abs(q);
         etemp := e;
         e := d;
         IF (abs(p) >= abs(0.5*q*etemp)) OR (p <= q*(a-x))
            OR (p >= q*(b-x)) THEN BEGIN
            IF x >= xm THEN
               e := a-x
            ELSE
               e := b-x;
            d := cgold*e
         END
         ELSE BEGIN
            d := p/q;
            u := x+d;
            IF (u-a < tol2) OR (b-u < tol2) THEN
               d := sign(tol1,xm-x)
         END
      END
      ELSE BEGIN
         IF x >= xm THEN
            e := a-x
         ELSE
            e := b-x;
         d := cgold*e
      END;
      IF abs(d) >= tol1 THEN
         u := x+d
      ELSE
         u := x+sign(tol1,d);
      fu := func(u);
      IF fu <= fx THEN BEGIN
         IF u >= x THEN a := x ELSE b := x;
         v := w;
         fv := fw;
         w := x;
         fw := fx;
         x := u;
         fx := fu
      END
      ELSE BEGIN
         IF u < x THEN a := u ELSE b := u;
         IF (fu <= fw) OR (w = x) THEN BEGIN
            v := w;
            fv := fw;
            w := u;
            fw := fu
         END
         ELSE IF (fu <= fv) OR (v = x) OR (v = w) THEN BEGIN
            v := u;
            fv := fu
         END
      END
   END;
   writeln('pause in routine BRENT - too many iterations');
99:
   xmin := x;
   brent := fx
END;

PROCEDURE mnbrak(VAR ax,bx,cx,fa,fb,fc: real);
{Pri zadanych AX, BX procedure hleda oblast minima funkce FUNC
 tim ze se pohybuje z kopce dolu doku nenarazi na protisvah.
 Pouziva parabolickou inter(extra)polaci k odhadu polohy minima.
 Na vystupu lezi BX mezi AX a CX a hodnota FUNC(BX) je nejmensi
 z danych 3 hodnot, ktere jsou  casti vystupu. Zvetseni intervalu
 se provadi nasobkem zlateho rezu GOLD. Parametr GLIMIT udava
 maximalni povolene zvetseni intervalu pri 1 kroku parabolicke
 extrpolace, TINY je minimalni hodnota jmenovatele pri parabolicke i(e).}
LABEL 99;
CONST
   gold = 1.618034;
   glimit = 100.0;
   tiny = 1.0e-20;
VAR
   ulim,u,r,q,fu,dum: real;

BEGIN
   fa := func(ax);
   fb := func(bx);
   IF fb > fa THEN BEGIN
      dum := ax;
      ax := bx;
      bx := dum;
      dum := fb;
      fb := fa;
      fa := dum
   END;
   cx := bx+gold*(bx-ax);
   fc := func(cx);
   WHILE fb >= fc DO BEGIN
      r := (bx-ax)*(fb-fc);
      q := (bx-cx)*(fb-fa);
      IF abs(q-r) > tiny THEN
         dum := abs(q-r)
      ELSE
         dum := tiny;
      IF q-r < 0.0 THEN dum := -dum;
      u := bx-((bx-cx)*q-(bx-ax)*r)/(2.0*dum);
      ulim := bx+glimit*(cx-bx);
      IF (bx-u)*(u-cx) > 0.0 THEN BEGIN
         fu := func(u);
         IF fu < fc THEN BEGIN
            ax := bx;
            fa := fb;
            bx := u;
            fb := fu;
            GOTO 99
         END
         ELSE IF fu > fb THEN BEGIN
            cx := u;
            fc := fu;
            GOTO 99
         END;
         u := cx+gold*(cx-bx);
         fu := func(u)
      END
      ELSE IF (cx-u)*(u-ulim) > 0.0 THEN BEGIN
         fu := func(u);
         IF fu < fc THEN BEGIN
            bx := cx;
            cx := u;
            u := cx+gold*(cx-bx);
            fb := fc;
            fc := fu;
            fu := func(u)
         END
      END
      ELSE IF (u-ulim)*(ulim-cx) >= 0.0 THEN BEGIN
         u := ulim;
         fu := func(u)
      END
      ELSE BEGIN
         u := cx+gold*(cx-bx);
         fu := func(u)
      END;
      ax := bx;
      bx := cx;
      cx := u;
      fa := fb;
      fb := fc;
      fc := fu
   END;
99:
END;

BEGIN
   LinminNcom := n;
   FOR j := 1 TO n DO BEGIN
      LinminPcom[j] := p[j];
      LinminXicom[j] := xi[j]
   END;
   ax := 0.0;
   xx := 1.0;
   mnbrak(ax,xx,bx,fa,fx,fb);
   fret := brent(ax,xx,bx,tol,xmin);
   FOR j := 1 TO n DO BEGIN
      xi[j] := xmin*xi[j];
      p[j] := p[j]+xi[j]
   END
END;

LABEL 99;
CONST
   itmax = 200;
   eps = 1.0e-10;
VAR
   j,its: integer;
   gg,gam,fp,dgg,sg: real;
   g,h,xi: ^RealArrayNP;
BEGIN
   new(g);
   new(h);
   new(xi);
   fp := fnc(p);
   dfnc(p,xi^);
   FOR j := 1 TO n DO BEGIN
      g^[j] := -xi^[j];
   END;
   Writeln(' g(0) - x y - ',g^[1],'  ',g^[2]);
   FOR its := 1 TO itmax DO BEGIN
      iter := its;
      linmin(p,g^,n,fret);
      IF 2.0*abs(fret-fp) <= ftol*(abs(fret)+abs(fp)+eps) THEN GOTO 99;
      fp := fnc(p);
      writeln('vysl iter',iter,'x y f',p[1],'  ',p[2],'  ',fp);
      dfnc(p,xi^);
      sg:=0.;
      FOR j := 1 TO n DO BEGIN
      g^[j] := -xi^[j];
      sg:=sg+g^[j]*g^[j]
      END;
      Writeln(' g(',iter:0,') - x y - ',g^[1],'  ',g^[2]);
      Readln;
      IF sg <= 0 THEN GOTO 99
   END;
   writeln('pause in routine FRPRMN');
   writeln('too many iterations');
   readln;
99:
   dispose(xi);
   dispose(h);
   dispose(g)
END;

VAR
    n,i,j,iter: Integer;
    fold, fnew, ftol : Real;
    pold,pnew :RealArrayNP;
LABEL 99;

BEGIN
    n:=2;
    WHILE true DO BEGIN
        WRITELN('ZADEJ X Y pocatecniho bodu a FTOL (|x|>0.9e30 konci)');
        READLN(pold[1],pold[2],ftol);
        IF abs(pold[1]) > 0.9E30 THEN GOTO 99;
        for i:=1 TO n DO pnew[i] := pold[i];
        fold := fnc(pold);
        writeln(' pocatek x y fold ',pnew[1],'  ',pnew[2],'  ',fold);
        frprmn(pnew,n,ftol,iter,fnew);
        writeln('iter =',iter);
        writeln(' x y fnew ',pnew[1],'  ',pnew[2],'  ',fnew);
   END;
99:
   Writeln(' Koncim... Dej Enter');
   Readln;
END.