function [lsg,wert,iteration,nfc,ngc,nGc]=trustregion(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
% ----------------------------------------------------------------------
% 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: 1000
% 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!
% nhc : 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 = [];
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 = 1000; 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);

while ~abbruch
   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 * 4; end;
   end;
   
   % Bestimme Suchrichtung s
   s = (H + nu*eye(n)) \ - g;
         
   % Bestimme den Grad der Uebereinstimmung zwischen f und der Approximation.
   xalt=x;
   hilfsf = feval(fktname,x+s);
   qs = f + g' * s + 0.5* s' * H * s;
   r = (f - hilfsf) / (f - qs);
   
   if r<0.25, nu = 4*nu; end;
   if r>0.75, nu = nu/2; end;
   if r>0,   x = x + s; end;
   
   falt = f; f = feval(fktname,x); normf = norm(f);
   galt = g; g = feval(gradname,x);normg = norm(g);
   iteration = iteration + 1;
   
   % Subroutine aufrufen fuer graphische Ausgabe, Minmax Bestimmung o.ae.
   if ~isempty(subroutine), feval(subroutine,x); end

   % 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;

lsg = x;
wert = f;
return;






