global numf numg numH geodata

geodata = struct('a',-1,'b',-1,'c',1,'d',1,'rho',2,'alpha',1,'beta',100);

% N = 1000;
N = 100;
% following is likely to have indefinite Hessian
z0 = ones(2*N,1);
% z0 = 2*rand(2*N,1)-1;
x = struct('p',z0);
nparams = struct('maxit',1000,'toler',1.0e-4,'method','chol');
sdparams = struct('maxit',1000,'toler',1.0e-4);

[inform,xnew] = SteepDescent(@geodesic,x,sdparams);
if inform.status == 0
  fprintf('CONVERGENCE FAILURE: %d steps were taken without\n', inform.iter);
  fprintf('gradient size decreasing below %8.4g.\n', sdparams.toler);
else
  fprintf('Success: %d steps taken\n', inform.iter);
end
[inform,xnew] = BFGS(@geodesic,x,nparams);
if inform.status == 0
  fprintf('CONVERGENCE FAILURE: %d steps were taken without\n', inform.iter);
  fprintf('gradient size decreasing below %8.4g.\n', nparams.toler);
else
  fprintf('Success: %d steps taken\n', inform.iter);
end
[inform,xnew] = Newton(@geodesic,x,nparams);
if inform.status == 0
  fprintf('CONVERGENCE FAILURE: %d steps were taken without\n', inform.iter);
  fprintf('gradient size decreasing below %8.4g.\n', nparams.toler);
else
  fprintf('Success: %d steps taken\n', inform.iter);
end
if N <= 5
  fprintf('  Ending point: '); fprintf('%8.4g ',xnew.p);
end
fprintf('\n  Ending value: '); fprintf('%8.4g\n',xnew.f);
fprintf('  No. function evaluations: %d',numf);
if N <= 5
  fprintf('\n  Ending gradient: '); fprintf('%8.4g ',xnew.g);
end
fprintf('\n  No. gradient evaluations: %d',numg);
fprintf('\n  Norm of ending gradient:  %8.4g\n\n\n', norm(xnew.g));

return;

[x,y] = meshgrid(geodata.a:0.1:geodata.c,geodata.b:0.1:geodata.d);
val = ones(size(x)) + ...
    geodata.alpha*exp(-geodata.beta*(x.^2 + y.^2));
contour(x,y,val);
figure(2); surf(x,y,val); figure(1);

hold on
m = length(xnew.p); n = m/2;
x = [geodata.a; xnew.p(1:2:m); geodata.c];
y = [geodata.b; xnew.p(2:2:m); geodata.d]; 
plot(x,y,'x-');
hold off