global numf numg numH geodata

x = struct('p',[-1.2; 1]);
nparams = struct('maxit',1000,'toler',1.0e-4,'method','chol');

[inform,xnew] = Newton(@obja,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
fprintf('  Ending point: '); fprintf('%8.4g ',xnew.p);
fprintf('\n  Ending value: '); fprintf('%8.4g ',xnew.f);
fprintf('; No. function evaluations: %d',numf);
fprintf('\n  Ending gradient: '); fprintf('%8.4g ',xnew.g);
fprintf('; No. gradient evaluations %d',numg);
fprintf('\n  Norm of ending gradient: %8.4g\n\n\n', norm(xnew.g));

[inform,xnew] = Newton(@objb,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
fprintf('  Ending point: '); fprintf('%8.4g ',xnew.p);
fprintf('\n  Ending value: '); fprintf('%8.4g ',xnew.f);
fprintf('; No. function evaluations: %d',numf);
fprintf('\n  Ending gradient: '); fprintf('%8.4g ',xnew.g);
fprintf('; No. gradient evaluations %d',numg);
fprintf('\n  Norm of ending gradient: %8.4g\n\n\n', norm(xnew.g));

[inform,xnew] = Newton(@objc,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
fprintf('  Ending point: '); fprintf('%8.4g ',xnew.p);
fprintf('\n  Ending value: '); fprintf('%8.4g ',xnew.f);
fprintf('; No. function evaluations: %d',numf);
fprintf('\n  Ending gradient: '); fprintf('%8.4g ',xnew.g);
fprintf('; No. gradient evaluations %d',numg);
fprintf('\n  Norm of ending gradient: %8.4g\n\n\n', norm(xnew.g));

x = struct('p',[1.1; 2.1]);
[inform,xnew] = Newton(@objd,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
fprintf('  Ending point: '); fprintf('%8.4g ',xnew.p);
fprintf('\n  Ending value: '); fprintf('%8.4g ',xnew.f);
fprintf('; No. function evaluations: %d',numf);
fprintf('\n  Ending gradient: '); fprintf('%8.4g ',xnew.g);
fprintf('; No. gradient evaluations %d',numg);
fprintf('\n  Norm of ending gradient: %8.4g\n\n\n', norm(xnew.g));

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

N = 1000;
% N = 100;
% N = 5;
% 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');
nparams = struct('maxit',1000,'toler',1.0e-6*sqrt(N),'method','direct');
sdparams = struct('maxit',1000,'toler',1.0e-6*sqrt(N));

fprintf('geodesic (%d) steepest descent\n', N);
[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
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));

fprintf('geodesic (%d) Newton\n', N);
[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));