function [lsg,wert,iteration,nfc,ngc,nGc]=newton(fktname,x,gradname,hessname,iterationmax,subroutine)
% ======================================================================
% Eine Implementation des Newton-Verfahrens mit Goldstein-Schrittlaenge
%
% Erstellt fuer die Veranstaltung
% Weiterbildenes Studium fuer Ingenieure WS 98/99
%
% S.F.
% ----------------------------------------------------------------------
% Aufrufparameter:
%  (lediglich die ersten vier Parameter muessen zwingend angegeben werden.)
% fktname : Name der Funktion F, deren Minimum bestimmt 
%   werden soll. An <fktname> wird ein Spaltenvektor(!) 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
% hessname : 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!
% ---------------------------------------------------------------------
% Programmvariablen:
% f : Wert der Funktion fname an der Stelle x.
% g : Wert des Gradienten von fname an der Stelle x.
% H : Wert der Hesse-Matrix von fname an der Stelle x.
% 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;
global nGc; nGc = 0;

verfahren='gold'; 	% Goldsteinschrittlaenge;

% fuer die eindimensionale Minimierung 
global fname; fname = fktname;
global xaktuell;
global gname; gname = gradname;
global s;



% Initialisierung
abbruch = 0;
iteration = 0;
lsg = [];
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 = 10; end;

% ======================================================================
% Das Newton-Verfahren
% ======================================================================
[g] = feval(gradname,x);
[f] = feval(fktname,x);
[H] = feval(hessname,x);

% 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 (Newton-Richtung)
  s = -H\g;

  % Schrittlaengenbestimmung 
  xaktuell=x;
  [alpha]=feval(verfahren,f,g);
  xalt=x;

  x = x + alpha * s;

  falt=f;  
  [f] = feval(fktname,x);
  normf = abs(f);
  galt=g;
  [g] = feval(gradname,x);
  normg = norm(g);
  [H] = feval(hessname,x);
  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 <= epsilon, abbruch = 1; end;
  if iteration >=iterationmax, abbruch = 1; end;

end;
% Rueckgabewerte festlegen - Ruecksprung zum aufrufenden Programm

lsg = x;
wert = f;
return;






