function [x, ithist, iflag] = shamnd( f, x, m, tolf, tolx, maxit )
%
%  function [x, ithist, iflag] = shamnd( f, x, m, tolf, tolx, maxit )
%
%  shamnd attempts to compute a root of F using Shamanskii's method
%
%  Input parameters:
%    f       name of a matlab function that evaluates 
%            f and its derivative.
%    x       initial iterate
%    m       number of iterations after which jacobian is recomputed
%            (m = 1: Newton's method; m = maxit : Chord method)
%    tolf    stopping tolerance (optional. Default tolf = 1.e-7)
%            Newton's method stops if  |f(x)| < tolf
%    tolx    stopping tolerance (optional. Default tolx = 1.e-7)
%            Newton's method stops if  |s| < tolx.
%    maxit   maximum number of iterations (optional. Default maxit = 100)
%
%
%  Output parameters:
%    x       approximation of the solution. 
%    ithist  array with the iteration history
%            The i-th row of ithist contains  [it, x, F, s]
%    ifag    return flag
%            iflag =  0  ||F(x)||_2 <= tolf 
%            iflag =  1  iteration terminated because maximum number of 
%                        iterations was reached. ||F(x)||_2 > tolf 
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  March 9, 2004
%

% set tolerances if necessary
if( nargin <= 4 ) tolf = 1.e-7; tolx = 1.e-7; maxit = 100; end
if( nargin <= 5 ) tolx = 1.e-7; maxit = 100; end
if( nargin <= 6 ) maxit = 100; end


   
it       = 0;
iflag    = 0;
norms    = 2*tolx;
[F,Jac]  = feval(f, x);
[L,U]    = lu(Jac);
   
while( it < maxit & norms > tolx & norm(F) > tolf )
   
     s     = - (U\(L\F));
     norms = norm(s);
   
     ithist(it+1,:) = [it, norm(x), norm(F), norms];
   
     x  = x+s;
     it = it+1;
     
     if( mod(it,m) == 0 )
         [F,Jac] = feval(f, x);
         [L,U]   = lu(Jac);
     else
         [F]     = feval(f, x);
     end
end


% check why the Shamanskii's method truncated and set iflag
if( norm(F) > tolf )
    % Shamanskii's method truncated because the maximum number of iterations
    % was reached
    iflag = 1;
    return
else
    % Shamanskii's method truncated because norm(F) <= tolf
    % print info for last iteration
    ithist(it+1,:) = [it, norm(x), norm(F),0];
end