## Copyright (C) 2013 Simone Panza
## Copyright (C) 2013 Fabio Cisaria
## Copyright (C) 2013 Carlo de Falco
##
## This file is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published
## by the Free Software Foundation; either version 3 of the License,
## or (at your option) any later version.
##
## This file is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License 
## along with this file; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {[@var{xi}, @var{x_iter}] =} @
## inextact_newton (@var{f}, @var{x0}, @var{tol}, @var{maxit}, @
## @var{eta}, @var{choice}, @var{iterative_method}, @var{restart})
##
## Solve the nonlinear equation @code{f(x) = 0} using the inexact
## Newton method, starting from the initial guess @var{x0}. 
##
## If an exact jacobian is not provided it is approximated by finite
## differences, the method for solving linear systems is specified in
## @var{iterative_method}, if not specified it defaults to 'gmres',
## other possible values are 'pcg' or 'bicgstab'. 
##
## Other input parameters have the following meanings
## @itemize @minus
## @item @var{f}: function handle to the function f(x)
## @item @var{x0}: initial guess
## @item @var{tol}: tolerance for termination test on residual 
## @item @var{maxit}: maximum number of nonlinear iterations
## @item @var{eta}: initial forcing term (must be in the interval 
## [0,1)), see [1] 
## @item @var{choice}: formula to use to select the forcing term may be
## 1 or 2, default is 1, see [1] for details.
## @item @var{iterative_method}: iterative method to solve the
## linearized system (default GMRES) 
## @item @var{restart}: restart parameter for the GMRES solver
## (ignored for other solvers)
## @end itemize
##
## The output parameters are:
## @itemize @minus
## @item @var{xi}: computed solution
## @iter @var{x_iter}: vector of the approximations of the solution at
## each step. 
## @end itemize
##	
## References:
## [1]  S.C. Eisenstat and H.F. Walker, "Choosing the Forcing
## Terms in an Inexact Newton Method." SIAM Journal on Scientific
## Computing, 17(1), pp. 16-16, 1996.
##
## [2] K. Turner and H. F. Walker, "Efficient high accuracy solutions
## with GMRES(m)." SIAM J. Sci. Statist. Comput.  13  (1992),  no. 3,
## 815--825. 
##
## @seealso{gmres, pcg, bicgstab, odebwenk}
## @end deftypefn


## Authors Simone Panza and Fabio Cisaria
## v1.1
## date 11-03-2013

function [x, x_iter] = inexact_newton (ff, x0, tol, maxit, eta, 
                                       choice, iterative_method,
                                       restart) 

  have_jacobian = false;
  if (iscell (ff) && numel (ff) == 2)
    f = ff{1};
    jac = ff{2};
    have_jacobian = true;
  elseif (is_function_handle (ff))
    f = ff;
  endif    
  
  if (nargin < 5 || isempty (eta))
    eta = .5; % default value, must be in the interval [0,1)
  endif
  
  eta = min (max (eta, 0), .99);
  
  if (nargin < 6 || isempty (choice))
    choice = 1; % default, see [1] for details
  endif
  
  if (nargin < 7 || isempty (iterative_method))
    iterative_method = 'gmres'; % default
    restart = [];
  endif

  if (nargin < 8)
    restart = [];
  endif
  
  switch iterative_method
    case 'gmres'
      solver = @(A, b, rtol) gmres (A, b, restart, rtol);
    case 'pcg'
      solver = @(A, b, rtol) pcg (A, b, rtol);
    case 'bicgstab'
      solver = @(A, b, rtol) bicgstab (A, b, rtol);
    otherwise
      error ('inexact_newton: invalid iterative solver %s', 
             iterative_method)  
  endswitch

  backtr_max = 10;    
  iter = 0;
  x = x_old = x_iter = x0;
  f_x_old = f_x = f(x);
  res = norm (f_x, 2);
  
  %% compute delta with the empirical formula given in [2]
  ord = 1; % 1st order
  delta = eps ^ (1 / (ord + 1)); 

  while ((res > tol) && (iter++ < maxit))

    if (! have_jacobian)     
      A = @(v) 1 / delta * (f (x + delta * v) - f_x);
    else
      A = jac (x);
    endif   

    x_old = x;
    f_x_old = f_x;
    
    [step, FLAG, RELRES, ITER, resvec] = solver (A, -f_x, eta); % use the chosen solver
    
    %% backtracking: reduce step size if nonlinear residual is bigger
    %% than the residual at the end of the past newton iteration 
    backtr_iter = 0;
    
    if 0
    factor = 1/2; % this is an arbitrary choice
    while ((norm (f(x + step), 2) > res) 
           && (backtr_iter++ < backtr_max))
      step *= factor;
      eta = eta/2 + 1/2; % see [1]
    endwhile
    
    if (backtr_iter == backtr_max)
      error('inexact_newton: residual cannot be reduced')
    endif
    end

    %% update x
    x = x + step;
    f_x = f (x);
    
    resold = res;
    res = norm (f_x, 2);
    x_iter = [x_iter x];
    
    %% compute the forcing term
    if (choice == 1) % see [1]
      
      %% eta = norm(F(x) - F(x_old) - F'(x_old)*step_old)/norm(F(x_old)) 
      %% F'(x_old)*step_old   can be approximated with first order
      %% finite differences formula. see [2] for details 
      
      eta_old = eta ^ ((1 + sqrt (5)) / 2);
      
      if (have_jacobian)
        eta = norm (f_x - f_x_old - A * step, 2) / resold;
      else
        eta = norm (f_x - f_x_old - A (step), 2) / resold;
      endif
      
    else % choice 2

      alpha = (1 + sqrt (5)) / 2; % see [1]
      gamma = 1;
      eta_old = gamma * eta ^ alpha;
      eta = gamma * (res / resold) ^ alpha;
      
    endif
      
    if (eta_old > .1) % safeguard
      eta = max (eta, eta_old);
    endif
    eta = min (eta, .99); % eta must belong to interval [0,1) 	
    
  endwhile
  
endfunction

%!test 
%! f = @(x) [((x(1)+1)^2 + 3*x(2)^7 + 7*(x(1)+1)); sum(x)+1];
%! sol = inexact_newton (f, [.1;.1], 1e-10, 100);
%! assert (sol, [-1; 0], 1e-9)
%! assert (f (sol), [0; 0], 1e-10)

%!test 
%! f = @(x) [((x(1)+1)^2 + 3*x(2)^7 + 7*(x(1)+1)); sum(x)+1];
%! jac = @(x) [2*(x(1)+1)+7, 3*7*x(2)^6; 1, 1];
%! [sol, it] = inexact_newton ({f, jac}, [.1;.1], 1e-10, 100);
%! assert (sol, [-1; 0], 1e-9)
%! assert (f (sol), [0; 0], 1e-10)

%!test 
%! f = @(x) [((x(1)+1)^2 + 3*x(2)^7 + 7*(x(1)+1)); sum(x)+1];
%! jac = @(x) [2*(x(1)+1)+7, 3*7*x(2)^6; 1, 1];
%! [sol, it] = inexact_newton ({f, jac}, [.1;.1], 1e-10, 100, .1, 2, 'bicgstab');
%! assert (sol, [-1; 0], 1e-9)
%! assert (f (sol), [0; 0], 1e-10)
