function [lsg,wert,iteration,nfc,ngc,nGc]=trustzei(fktname,x,gradname,hessename,iterationmax,subroutine)
% ======================================================================
% Eine Implementation eines Trust-Region Verfahrens
%
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% Helge Baumann, 09.10.1998
% Aenderungen S.F.
% ----------------------------------------------------------------------
% Aurufparameter:
%  (lediglich die ersten vier Parameter muessen zwingend angegeben werden.)
% fktname : Name der Funktion F, deren Minimum bestimmt 
%   werden soll. An <fktname> wird ein Spaltenwektor(!) uebergeben. 
%   Es wird eine reelle Zahl als Bild erwartet.
% x : aktueller Iterationspunkt. 
% gradname : Name der Funktion G, die den Gradienten von F als 
%   Spaltenvektor zurueckgibt
% hessename : Name der Funktion H, die die Hessematrix von F zurueckgibt.
% iterationmax : Anzahl der maximalen Iterationen. Standard: 5
% subroutine : Name einer Unterroutine, die in jedem Iterationsschritt
%   aufgerufen wird. Dient zum Protokolieren des Iterationsprozesses.
%   Wird <''> uebergeben, erfolgt kein Aufruf. Standard: ''
% ---------------------------------------------------------------------
% Rueckgabeparameter:
% lsg : Loesung (Spaltenvektor!)
% wert: Zielfunktionswert an der Stelle lsg
% iteration : Anzahl der benoetigten Iterationen
% nfc : Anzahl der Funktionsauswertungen, muss allerdings in der
%   Unterroutine <fktname> gesetzt werden!
% ngc : Anzahl der Gradientenauswertungen, muss allerdings in der
%   Unterroutine <gradname> gesetzt werden!
% nGc : Anzahl der Hessematrizenauswertungen, muss allerdings in der
%   Unterroutine <hessename> gesetzt werden!
% ---------------------------------------------------------------------
% Programmvariablen:
% f : Wert der Funktion fname an der Stelle x.
% g : Wert des Gradienten von fname an der Stelle x.
% ======================================================================
global nfc; nfc = 0;
global ngc; ngc = 0;
global nGc; nGc = 0;
 
% Initialisierung
epsilon = 1E-10; 
abbruch = 0;
iteration = 0;
lsg = [];

% Parameter:	nach Algorithmus 8.4
%	0 < o1 < 1 < o2
o1=0.5;
o2=4;
%	0 <= rho0 < rho1 < rho2 < 1
rho0=0;
rho1=0.25;
rho2=0.75;

nu = 1;
[n,m] = size(x);

% Konstanten fuer die Abbruchbedingung
tf = 10^(-8); tf2 = sqrt(tf); tf3 = tf^(1/3);

if nargin < 5, subroutine = ''; end;
if nargin < 4, iterationmax = 5; end;
 
% Subroutine aufrufen fuer graphische Ausgabe, Minmax Bestimmung o.ae.
if ~isempty(subroutine), feval(subroutine,x); end
 
% ======================================================================
% Der Trust-Region-Algorithmus
% ======================================================================
[f] = feval(fktname,x);
[g] = feval(gradname,x);

%	Algerithmus 8.4

while ~abbruch	

%  %	While g'(x)+nu*I nicht positiv definit ...

   H = feval(hessename,x);
   p = 1;
   
   while (p>0)
      % Test auf positive Definitheit mittels Cholesky Zerlegung.
      [C,p] = chol(H + nu * eye(n));
      % Wenn p>0 ist die Matrix nicht positiv definit.
      if p>0, nu = nu * o2; end;
   end;
   
%  % Bestimme Suchrichtung s (aus 8.4)
   s = (H + nu*eye(n)) \ - g;
         
%  % Bestimme den Grad der Uebereinstimmung zwischen f und der Approximation.
   % (aus 8.3)
   xalt=x;
   hilfsf = feval(fktname,x+s);
   qs = f + g' * s + 0.5* s' * H * s;
   r = (f - hilfsf) / (f - qs);

   nualt=nu;	%	Speichern von nu fuer Zeichnung

%  %	Berechen von neuem nu und neuem x
   if r<rho1, nu = o2*nu; end;
   if r>rho2, nu = o1*nu; end;
   if r>rho0,  x = x + s; end;

%  %    Neue Funktionswert, Wert des Gradienten, Abbruchbedingung

   falt = f; f = feval(fktname,x); normf = norm(f);
   galt = g; g = feval(gradname,x);normg = norm(g);
   iteration = iteration + 1;

%-------------------- Zeichnung, gehoert nicht zum Algorithmus --------%
   if r>0, 	%	Zeichnen                                       %
	figure                                                         %
	h=norm(s); 	%	Trustregieon-Radius                    %
	links=xalt(1)-h;rechts=xalt(1)+h;                              %
	oben=xalt(2)+h;unten=xalt(2)-h;                                %
	axis([links,rechts,unten,oben])                                %
	hold on                                                        %
	xpoints=61;ypoints=61;                                         %
	xrange=linspace(links,rechts,xpoints);                         %
	yrange=linspace(unten,oben,ypoints);                           %
 	[xgitter,ygitter] = meshgrid(xrange,yrange);                   %
	for i=1:xpoints                                                %
      	  [zgitter(i,1:ypoints)] = ...                                 %
	     feval(fktname,[xrange(i)*ones(1,ypoints);yrange(1:ypoints)]);
	end                                                            %
 	% Erzeugzung der Hoehenlinien   VON F                          %
	flines=[linspace(f,falt,5),linspace(falt,10*falt,5)];          %
   	contour(xgitter,ygitter,zgitter',flines);                      %
        % Hoehenlienen quadratische Hilfsfunktion                      %
	for i=1:xpoints                                                %
	   for j=1:ypoints                                             %
		st=([xgitter(i,j);ygitter(i,j)]-xalt);                 %
   		qgitter(i,j) = falt + galt' * st + 0.5* st' * H * st;  %      
      	   end;                                                        %
	end;                                                           %
   	contour(xgitter,ygitter,qgitter,flines,'k-.');                 %
	xlabel(['nu=',num2str(nualt)])                                 %
	plot(xalt(1),xalt(2),'r*');                                    %
	plot(x(1),x(2),'g+')                                           %
	 hold off                                                      %
   end;                                                                %
   % Subroutine aufrufen fuer graphische Ausgabe, Minmax Bestimmung o.ae.
   if ~isempty(subroutine), feval(subroutine,x); end                   %
                                                                       %
% ----------------------Ende Zeichnung --------------------------------%

   % Abbruchbedingungen ueberpruefen
   if ((falt - f) <= (tf * (1+normf))) & ...
      ( norm(xalt-x) <= (tf2 * (1 + norm(x)))) & ...
      ( normg <= (tf3 * (1+normf)) ) 
     abbruch = 1; 
   end;
   if normg <= realmin, abbruch = 1; end;
   if iteration >=iterationmax, abbruch = 1; end;
end;	%	ENDE WHILE

lsg = x;
wert = f;
return;












