{ukazka Powellovy metody konjugovanych smeru,
konverguje rychle, zde extremne rychle nebot se hleda minimum
kvadraticke formy
}
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 ZadSmer(n: integer; VAR xi: RealArrayNPbyNP);
VAR
   i,j : Integer;
BEGIN
     FOR i:=1 TO n DO BEGIN
            FOR j:=1 TO n DO  xi[i,j] := 0;
            xi[i,i] :=1
            END;
END;


PROCEDURE powell(VAR p: RealArrayNP;
                VAR xi: RealArrayNPbyNP;
                     n: integer;
                  ftol: real;
              VAR iter: integer;
              VAR fret: real);
{Hleda minimum funkce FNC(X) N promennych (X je vektor dimenze N)
 Pocatecni bod je P[1..N], XI[1..N,1..N] je matice pocatecni smeru,
 obvykle N jednotkovych vektoru a FTOL je pomerna tolerance funkcni
 hodnoty minima. Neschpnost zmensit hodnotu funkce o vice nez tuto
 hodnotu znamena hotovo. Na vystupu P je nejlepsi nalezeny bod,
 XI je soucasna sada smeru, FRET je hodnota funkce je bode P,
 ITER je pocet provedenych iteraci. Procedura pouziva LINMIN
 pro hledani minima ve smeru.}

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;
VAR
   j,ibig,i: integer;
   t,fptt,fp,del: real;
   pt,ptt,xit: ^RealArrayNP;
BEGIN
   new(pt);
   new(ptt);
   new(xit);
   fret := fnc(p);
   FOR j := 1 TO n DO pt^[j] := p[j];
   iter := 0;
   WHILE true DO BEGIN
      iter := iter+1;
      writeln('probiha iterace  ',iter);
      fp := fret;
      ibig := 0;
      del := 0.0;
      FOR i := 1 TO n DO BEGIN
         FOR j := 1 TO n DO
            xit^[j] := xi[j,i];
         fptt := fret;
         writeln('smer ',i,'sx sy',xit^[1],'  ',xit^[2]);
         linmin(p,xit^,n,fret);
         writeln('vysl i',i,'x y f',p[1],'  ',p[2],'  ',fret);
         IF abs(fptt-fret) > del THEN BEGIN
            del := abs(fptt-fret);
            ibig := i
         END
      END;
      IF 2.0*abs(fp-fret)
         <= ftol*(abs(fp)+abs(fret)) THEN GOTO 99;
      IF iter = itmax THEN BEGIN
         writeln('pause in routine POWELL');
         writeln('too many interations');
         readln
      END;
      FOR j := 1 TO n DO BEGIN
         ptt^[j] := 2.0*p[j]-pt^[j];
         xit^[j] := p[j]-pt^[j];
         pt^[j] := p[j]
      END;
      fptt := fnc(ptt^);
      IF fptt < fp THEN BEGIN
         t := 2.0*(fp-2.0*fret+fptt)*sqr(fp-fret-del)-del*sqr(fp-fptt);
         IF t < 0.0 THEN BEGIN
            writeln('delam 2 minimalizaci','sx sy',xit^[1],'  ',xit^[2]);
            linmin(p,xit^,n,fret);
            writeln('vysledek x y f',p[1],'  ',p[2],'  ',fret);
            FOR j := 1 TO n DO
               xi[j,ibig] := xit^[j]
         END
      END;
   Readln;
   END;
99:
   dispose(xit);
   dispose(ptt);
   dispose(pt)
END;

VAR
    n,i,j,iter: Integer;
    fold, fnew, ftol : Real;
    pold,pnew :RealArrayNP;
    xold :RealArrayNPbyNP;
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);
        zadsmer(n,xold);
        powell(pnew,xold,n,ftol,iter,fnew);
        writeln('iter =',iter);
        writeln(' x y fnew ',pnew[1],'  ',pnew[2],'  ',fnew);
   END;
99:
   Writeln(' Koncim... Dej Enter');
   Readln;
END.