function rosenbrock_dogleg

% set to 1 for a Cauchy-point based algorithm; otherwise dogleg is used
use_Cauchy = 0;

x0=[-1.2; 1];
xs=[1;1];


% Trust-region/dogleg algorithm
Delta = 1.0;
xk    = x0;
df0   = df(x0);

Nmax  = 1000;

for i = 1:Nmax,
    g =  df(xk);
    B = df2(xk);
    %
    printf('%3d: f=%15.6e, |df|=%15.6e, |x-x*| = %15.6e Delta = %15.6e\n', ...
           i, f(xk), norm(g), norm(xk-xs), Delta);
    %
    if norm(g) <= 1.0E-12*norm(df0)
        % converged!
        break;
    end
    % compute Cauchy point
    pC = cauchy(g,B,Delta);
    if min(eig(B)) <= 0 || use_Cauchy
        % B not positive definite - use Cauchy point step
        p = pC;
    else
        % else - dogleg
        p = dogleg(g,B,Delta,pC);
    end
    %
    %  check the actual reduction
    rho = (f(xk) - f(xk + p)) / (-g'*p - p'*B*p/2);
    % update the trust-region radius
    if rho < 0.25
        Delta = 0.25*Delta;
    elseif rho > 0.75 && Delta-norm(p) >= -1.0E-06*Delta
        Delta = min(2*Delta, 2.0);
    end
    % accept the step, if the ratio is OK
    if rho > 1.0E-03
        xk = xk+p;
    end
end

function pC = cauchy(g,B,Delta)
if g'*B*g <= 0
    pC = -Delta/norm(g)*g;
else
    alphaS = g'*g/(g'*B*g);
    if norm(alphaS*g) <= Delta
        pC = -alphaS*g;
    else
        pC = -Delta/norm(g)*g;
    end
end

function p = dogleg(g,B,Delta,pC)
if norm(pC) >= Delta
    p = pC;
else
    pB = -B\g;
    if norm(pB) <= Delta
        p = pB;
    else
        % solve the eqn |pC + alpha*(pB-pC)|^2-Delta^2 = 0
        a  = norm(pB-pC)^2;
        b  = pC'*(pB-pC);
        c  = norm(pC)^2 - Delta^2;
        D  = b^2 - a*c;
        alpha = (-b + sqrt(D))/a;
        p  = pC + alpha*(pB-pC);
    end
end

function r=f(x)
r = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

function r=df(x)
r=[-400*(x(2)-x(1)^2)*x(1) - 2*(1-x(1)); 200*(x(2)-x(1)^2)];
    
function r=df2(x)
r=[[400*(3*x(1)^2-x(2))+2, -400*x(1)];  [-400*x(1), 200]];