function [lsg,wert,iteration,nfc,ngc]=bfgs(fktname,x,gradname,iterationmax,subroutine)
% ======================================================================
% Eine Implementation des BFGS-Algorithmus mit der 
% Schrittlaengenbestimmung nach Goldstein
%
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% Helge Baumann, 10.9.1998
% ----------------------------------------------------------------------
% Aurufparameter:
%  (lediglich die ersten drei 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
% 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!
% ---------------------------------------------------------------------
% Programmvariablen:
% f : Wert der Funktion fname an der Stelle x.
% g : Wert des Gradienten von fname an der Stelle x.
% H : Quasi-Newton-Matrix
% normf, normg : Die Norm von f bzw. g.
% xalt,falt,galt : Die Werte von x bzw. f,g im vorherigen Schritt.
% alpha : Schrittlaenge
% s : Suchrichtung
% p,y,pty,Hy : Hilfsvektoren bzw. Hilfsvariable
% iterationmax : Anzahl der maximalen Iterationen
% abbruch, tf,tf2,tf3  :  Hilfsvariable
% ======================================================================
global nfc; nfc = 0;
global ngc; ngc = 0;
 
 % fuer die Eindimensionale Minimierung in der Funktion phi:
 global fname; fname = fktname;
 global xaktuell
 global s
 
% Initialisierung
abbruch = 0;
iteration = 0;
lsg = [];
rho = 0.1;
epsilon = 10^(-8);
format long

% 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;
 
% ======================================================================
% Das BFGS-Verfahren
% ======================================================================
[g] = feval(gradname,x);
[f] = feval(fktname,x);
H = eye(size(x,1));
% Subroutine aufrufen fuer graphische Ausgabe, Minmax Bestimmung o.ae.
if ~isempty(subroutine), feval(subroutine,x); end

% Ueberpruefung, ob man bereits in einem stationaeren Punkt startet
if norm(g) <= epsilon, abbruch = 1; end;

while (abbruch ~= 1)

  % Suchrichtung 
  s = - H * g;
  
  % Bestimmung der Schrittlaenge alpha:
  % Schrittlaengenstrategie nach Goldstein:
  alpha = 0; beta=1;
  xaktuell = x;  % wird fuer phi benoetigt.
  phi0 = phi(0); phistrich0 = g' * s;

  while (phi(beta) < phi0+(1-rho)*beta*phistrich0 )
    alpha = beta; beta = 2*beta;
  end;
  if phi(beta) > phi0 + rho*beta*phistrich0
    gamma = 0;           % damit in nachfolgender Schleife nicht vor dem 
    phigamma =  phi0+1;   % ersten Durchlauf abgebrochen wird.
    while ( (phigamma < phi0 + (1-rho)*gamma*phistrich0) | ...
            (phigamma > phi0 + rho*gamma*phistrich0 ))
	gamma = (alpha+beta) / 2;
	phigamma = phi(gamma);
	if phigamma < phi0 + (1-rho)*gamma * phistrich0, alpha=gamma; end;
	if phigamma > phi0 + rho * gamma * phistrich0, beta = gamma; end;
    end;     % end of while
    alpha = gamma;
  else
    alpha = beta;
  end;
  % in alpha steht nun eine Schrittlaenge, 
  % die den Goldsteinbedingungen gennuegt.
 
  xalt = x;
  x = x + alpha * s;
  
  falt = f;
  [f] = feval(fktname,x);
  normf = abs(f);
  galt = g;
  [g] = feval(gradname,x);
  normg = norm(g);
  iteration = iteration +1;

  % Aufdatierung von H nach BFGS-Formel:
  p = x - xalt;
  y = g-galt;
  pty = p' * y;
  % Notbremse falls Nenner Null wird.
  if abs(pty) < eps
     error('Abbruch -- Nenner in Aufdatierungsformel ist Null !!!');
  end;
  Hy = H*y;
  H = H + (1+ y'*Hy/pty) * p*p'/pty - (p*Hy'+Hy*p')/pty;

  % 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 <= epsilon, abbruch = 1; end;
  if iteration >=iterationmax, abbruch = 1; end;

end;

% Rueckgabewerte festlegen - Ruecksprung zum aufrufenden Programm

lsg = x;
wert = f;
return;






