function [xmin,fmin] = Newton_backtrack(x,fname,prtout,tol)
% Call : [xmin,fmin] = Newton_backtrack(x,fname,prtout,tol)
%
% Written by Emre Mengi - April 17, 2008
% Last Modified - May 29, 2008
%
% Task :
% Applies Newton Iteration with backtracking to find a root of
% the function fname starting from x.
%
% Input parameters:
% x	(vector)	   : Initial point for Newton's method.
% fname (string)	   : Name of the function for which a root is
%			     sought. This function is supposed to take one input
%			     argument, the point where the function should be
%			     evaluated. It returns two output arguments, the
%			     function value (possibly a vector) and the Jacobian.
% prtout (integer)	   : If prtout > 0 (Default), print the iterate and the function
%			     value at each iteration. 
%			     If prtout <= 0, don't print.
% tol (real scalar)	   : Tolerance for termination. When the 2-norm of the function
%			     value is at most tol, stop. (Default --- tol=10^-6)
%
% Output arguments:
% xmin (vector)		   : Computed root.
% fmin (vector)		   : Function value at xmin; cannot be greater than
%			     tol in 2-norm. For better accuracy set
%			     tol small. Setting tol too small causes convergence
%			     problems due to rounding errors.


if nargin < 2
	error('two input arguments - initial point and fun name - required')
end

if nargin < 3
	prtout = 1; % default print
end

if nargin < 4
	tol = 10^-6;
end

mu = 1/4;
iter = 0;
[fval, J] = feval(fname,x);
if prtout > 0
	fprintf('iter %3d\n', iter);
	fprintf('x =');
	fprintf(' %.9f', x);
	fprintf('\nf =');
	fprintf(' %.9f', fval);
	fprintf('\n');
end

while abs(norm(fval)) > tol & iter < 500
	
	p = -J\fval;
	
	innerit = 0;
	normfval = norm(fval);
	alpha = 1;	
	xt = x + alpha*p;
	[fvaln, Jn] = feval(fname,xt);
	  
	while (innerit < 20) & (norm(fvaln) > (1 - alpha*mu)*normfval) 
		alpha = alpha/2;
		xt = x + alpha*p;
		[fvaln,Jn] = feval(fname,xt);
		innerit = innerit + 1; 
	end	
	
	x = xt;
	fval = fvaln;
	J = Jn;
	
	iter = iter + 1;
	
	if prtout > 0
		fprintf('iter %3d\n', iter);
		fprintf('x =');
		fprintf(' %.9f', x);
		fprintf('\nnormf =');
		fprintf(' %.9f', normfval);
		fprintf('\n');
	end

end

xmin = x;
fmin = fval;

if prtout > 0
	fprintf('\nNumber of Newton iterations : %d\n',iter);
	fprintf('Computed root : \n');
	fprintf(' %.10f', xmin);
	fprintf('\nfunction value at termination : \n');
	fprintf(' %.10f', fmin);
	fprintf('\nnorm of function at termination : \n');
	fprintf(' %.10f', norm(fmin));
	fprintf('\n');
end

if iter >= 500
	fprintf('\n');
	warning('Newtons method could not converge after 500 iterations') 
end

return