## Copyright (C) 2009-2011 Carlo de Falco
##  
## This file is part of Octave.
## 
## Octave 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.
## 
## Octave 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 Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
## @deftypefnx {Function File} {@var{x} =} gmres (@var{A}, @var{b}, @var{m}, @var{rtol}, @var{maxit}, @var{P})
## @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} gmres (@dots{})
## Solve @code{A x = b} using the Preconditioned GMRES iterative method
## with restart, a.k.a. PGMRES(m).
##
## @itemize @minus
## @item @var{rtol} is the relative tolerance,
## if not given or set to [] the default value 1e-6 is used.
##
## @item @var{maxit} is the maximum number of outer iterations,
## if not given or set to [] the default value
## @code{min (10, numel (b) / restart)} is used.
##
## @item @var{x0} is the initial guess,
## if not given or set to [] the default value @code{zeros(size (b))} is used.
##
## @item @var{m} is the restart parameter,
## if not given or set to [] the default value @code{numel (b)} is used.
## @end itemize
##
## Argument @var{A} can be passed as a matrix, function handle, or
## inline function @code{f} such that @code{f(x) = A*x}.
##
## The preconditioner @var{P} is given as @code{P = M1 * M2}.
## Both @var{M1} and @var{M2} can be passed as a matrix, function handle, or
## inline function @code{g} such that @code{g(x) = M1\x} or @code{g(x) = M2\x}.
##
## Besides the vector @var{x}, additional outputs are:
##
## @itemize @minus
## @item @var{flag} indicates the exit status:
##
## @table @asis
##   @item 0 : iteration converged to within the specified tolerance
##
##   @item 1 : maximum number of iterations exceeded
##
##   @item 2 : unused, but skipped for compatibility
##
##   @item 3 : algorithm reached stagnation
## @end table
##
## @item @var{relres} is the final value of the relative residual.
##
## @item @var{iter} is a vector containing the number of outer iterations and
## total iterations performed.
##
## @item @var{resvec} is a vector containing the relative residual at each
## iteration.
## @end itemize
##
## @seealso{pcg, cgs, bicgstab}
## @end deftypefn

function [x, flag, prec_res_norm, itcnt] = gmres (A, b, restart, rtol, maxit, M1, M2, x0)

  if (nargin < 2 || nargin > 8)
    print_usage ();    
  end
      
  if (ischar (A))
    Ax = str2func (A);
  elseif (ismatrix (A))
    Ax = @(x) A*x;
  elseif (isa (A, "function_handle"))
    Ax = A;
  else
    error ("gmres: A must be a function or matrix");
  endif

  if (nargin < 3 || isempty (restart))
    restart = rows (b);
  endif

  if (nargin < 4 || isempty (rtol))
    rtol = 1e-6;
  endif

  if (nargin < 5 || isempty (maxit))
    maxit = min (rows (b)/restart, 10);
  endif

  if (nargin < 6 || isempty (M1))
    M1m1x = @(x) x;
  elseif (ischar (M1))
    M1m1x = str2func (M1);
  elseif (ismatrix (M1))
    M1m1x = @(x) M1 \ x;
  elseif (isa (M1, "function_handle"))
    M1m1x = M1;
  else
    error ("gmres: preconditioner M1 must be a function or matrix");
  endif
  
  if (nargin < 7 || isempty (M2))
    M2m1x = @(x) x;
  elseif (ischar (M2))
    M2m1x = str2func (M2);
  elseif (ismatrix (M2))
    M2m1x = @(x) M2 \ x;
  elseif (isa (M2, "function_handle"))
    M2m1x = M2;
  else
    error ("gmres: preconditioner M2 must be a function or matrix");
  endif

  Pm1x = @(x) M2m1x (M1m1x (x));

  if (nargin < 8 || isempty (x0))
    x0 = zeros (size (b));
  endif

  x_old = x0; 
  x = x_old;
  prec_res = Pm1x (b - Ax (x_old));
  prec_res_norm = norm (prec_res, 2);
  
  B = zeros (restart + 1, 1);
  V = zeros (rows (x), restart);
  H = zeros (restart + 1, restart);

  ## begin loop
  iter = 1;
  restart_it  = restart + 1; 
  resids      = zeros (maxit, 1);
  resids(1)   = prec_res_norm;
  prec_b_norm = norm (Pm1x (b), 2);
  flag        = 1;

  while ((iter <= maxit * restart) && (prec_res_norm > rtol * prec_b_norm))
  
    ## restart
    if (restart_it > restart)
      restart_it = 1;
      x_old = x;	      
      prec_res = Pm1x (b - Ax (x_old));
      prec_res_norm = norm (prec_res, 2);
      B(1) = prec_res_norm;
      H(:) = 0;
      V(:, 1) = prec_res / prec_res_norm;
    endif  
    
    ## basic iteration
    tmp = Pm1x (Ax (V(:, restart_it)));
    [V(:,restart_it+1), H(1:restart_it+1, restart_it)] = mgorth (tmp, V(:,1:restart_it));
    
    Y = (H(1:restart_it+1, 1:restart_it) \ B (1:restart_it+1));
	      
    little_res = B(1:restart_it+1) - H(1:restart_it+1, 1:restart_it) * Y(1:restart_it);
    prec_res_norm = norm (little_res, 2);
	      
    x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
    
    resids(iter) = prec_res_norm ;
    if (norm (x - x_old, inf) <= eps)
      flag = 3;
      break
    endif

    restart_it++ ; iter++;
  endwhile

  if (prec_res_norm > rtol * prec_b_norm)
    flag = 0;
  endif

  resids = resids(1:iter-1);
  itcnt = [floor(maxit/restart), rem(maxit, restart)];
endfunction


%!shared A, b, dim
%! dim = 100;
%!test
%! A = spdiags ([-ones(dim,1) 2*ones(dim,1) ones(dim,1)], [-1:1], dim, dim);
%! b = ones(dim, 1); 
%! x = gmres (A, b, 10, 1e-10, dim, @(x) x./diag(A), [],  b);
%! assert(x, A\b, 1e-9*norm(x,inf));
%!
%!test
%! x = gmres (A, b, dim, 1e-10, 1e4, @(x) diag(diag(A))\x, [],  b);  
%! assert(x, A\b, 1e-7*norm(x,inf));
%!
%!test
%! A = spdiags ([[1./(2:2:2*(dim-1)) 0]; 1./(1:2:2*dim-1); [0 1./(2:2:2*(dim-1))]]', -1:1, dim, dim);
%! A = A'*A;
%! b = rand (dim, 1);
%! [x, resids] = gmres (@(x) A*x, b, dim, 1e-10, dim, @(x) x./diag(A), [],  []);
%! assert(x, A\b, 1e-9*norm(x,inf))
%! x = gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) diag(diag(A))\x, [],  []);
%! assert(x, A\b, 1e-9*norm(x,inf));
%!test
%! x =  gmres (@(x) A*x, b, dim, 1e-10, 1e6, @(x) x./diag(A), [],  []);
%! assert(x, A\b, 1e-7*norm(x,inf));