% -*- texinfo -*-
% @deftypefn {Function File} {@var{x} =} pcg (@var{a}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
% @deftypefnx {Function File} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@dots{})
%
% Solves the linear system of equations @code{@var{a} * @var{x} =
% @var{b}} by means of the Preconditioned Conjugate Gradient iterative
% method. The input arguments are
%
% @itemize
% @item
% @var{a} can be either a square (preferably sparse) matrix or a
% function handle, inline function or string containing the name
% of a function which computes @code{@var{a} * @var{x}}. In principle
% @var{a} should be symmetric and positive definite; if @code{pcg}
% finds @var{a} to not be positive definite, you will get a warning
% message and the @var{flag} output parameter will be set.
%
% @item
% @var{b} is the right hand side vector.
%
% @item
% @var{tol} is the required relative tolerance for the residual error,
% @code{@var{b} - @var{a} * @var{x}}. The iteration stops if @code{norm
% (@var{b} - @var{a} * @var{x}) <= @var{tol} * norm (@var{b} - @var{a} *
% @var{x0})}. If @var{tol} is empty or is omitted, the function sets
% @code{@var{tol} = 1e-6} by default.
%
% @item
% @var{maxit} is the maximum allowable number of iterations; if
% @code{[]} is supplied for @code{maxit}, or @code{pcg} has less
% arguments, a default value equal to 20 is used.
%
% @item
% @var{m} = @var{m1} * @var{m2} is the (left) preconditioning matrix, so that the iteration is
% (theoretically) equivalent to solving by @code{pcg} @code{@var{P} *
% @var{x} = @var{m} \ @var{b}}, with @code{@var{P} = @var{m} \ @var{a}}.
% Note that a proper choice of the preconditioner may dramatically
% improve the overall performance of the method. Instead of matrices
% @var{m1} and @var{m2}, the user may pass two functions which return
% the results of applying the inverse of @var{m1} and @var{m2} to
% a vector (usually this is the preferred way of using the preconditioner).
% If @code{[]} is supplied for @var{m1}, or @var{m1} is omitted, no
% preconditioning is applied. If @var{m2} is omitted, @var{m} = @var{m1}
% will be used as preconditioner.
%
% @item
% @var{x0} is the initial guess. If @var{x0} is empty or omitted, the
% function sets @var{x0} to a zero vector by default.
% @end itemize
%
% The arguments which follow @var{x0} are treated as parameters, and
% passed in a proper way to any of the functions (@var{a} or @var{m})
% which are passed to @code{pcg}. See the examples below for further
% details. The output arguments are
%
% @itemize
% @item
% @var{x} is the computed approximation to the solution of
% @code{@var{a} * @var{x} = @var{b}}.
%
% @item
% @var{flag} reports on the convergence. @code{@var{flag} = 0} means
% the solution converged and the tolerance criterion given by @var{tol}
% is satisfied. @code{@var{flag} = 1} means that the @var{maxit} limit
% for the iteration count was reached. @code{@var{flag} = 3} reports that
% the (preconditioned) matrix was found not positive definite.
%
% @item
% @var{relres} is the ratio of the final residual to its initial value,
% measured in the Euclidean norm.
%
% @item
% @var{iter} is the actual number of iterations performed.
%
% @item
% @var{resvec} describes the convergence history of the method.
% @code{@var{resvec} (i,1)} is the Euclidean norm of the residual, and
% @code{@var{resvec} (i,2)} is the preconditioned residual norm,
% after the (@var{i}-1)-th iteration, @code{@var{i} =
% 1, 2, @dots{}, @var{iter}+1}. The preconditioned residual norm
% is defined as
% @code{norm (@var{r}) ^ 2 = @var{r}' * (@var{m} \ @var{r})} where
% @code{@var{r} = @var{b} - @var{a} * @var{x}}, see also the
% description of @var{m}. If @var{eigest} is not required, only
% @code{@var{resvec} (:,1)} is returned.
%
% @item
% @var{eigest} returns the estimate for the smallest @code{@var{eigest}
% (1)} and largest @code{@var{eigest} (2)} eigenvalues of the
% preconditioned matrix @code{@var{P} = @var{m} \ @var{a}}. In
% particular, if no preconditioning is used, the estimates for the
% extreme eigenvalues of @var{a} are returned. @code{@var{eigest} (1)}
% is an overestimate and @code{@var{eigest} (2)} is an underestimate,
% so that @code{@var{eigest} (2) / @var{eigest} (1)} is a lower bound
% for @code{cond (@var{P}, 2)}, which nevertheless in the limit should
% theoretically be equal to the actual value of the condition number.
% The method which computes @var{eigest} works only for symmetric positive
% definite @var{a} and @var{m}, and the user is responsible for
% verifying this assumption.
% @end itemize
%
% Let us consider a trivial problem with a diagonal matrix (we exploit the
% sparsity of A)
%
% @example
% @group
% n = 10;
% a = diag (sparse (1:n));
% b = rand (n, 1);
% [l, u, p, q] = luinc (a, 1.e-3);
% @end group
% @end example
%
% @sc{Example 1:} Simplest use of @code{pcg}
%
% @example
% x = pcg(A,b)
% @end example
%
% @sc{Example 2:} @code{pcg} with a function which computes
% @code{@var{a} * @var{x}}
%
% @example
% @group
% function y = apply_a (x)
% y = [1:N]'.*x;
%
%
% x = pcg ('apply_a', b)
% @end group
% @end example
%
% @sc{Example 3:} @code{pcg} with a preconditioner: @var{l} * @var{u}
%
% @example
% x = pcg (a, b, 1.e-6, 500, l*u);
% @end example
%
% @sc{Example 4:} @code{pcg} with a preconditioner: @var{l} * @var{u}.
% Faster than @sc{Example 3} since lower and upper triangular matrices
% are easier to invert
%
% @example
% x = pcg (a, b, 1.e-6, 500, l, u);
% @end example
%
% @sc{Example 5:} Preconditioned iteration, with full diagnostics. The
% preconditioner (quite strange, because even the original matrix
% @var{a} is trivial) is defined as a function
%
% @example
% @group
% function y = apply_m (x)
% k = floor (length (x) - 2);
% y = x;
% y(1:k) = x(1:k)./[1:k]';
%
%
% [x, flag, relres, iter, resvec, eigest] = ...
% pcg (a, b, [], [], 'apply_m');
% semilogy (1:iter+1, resvec);
% @end group
% @end example
%
% @sc{Example 6:} Finally, a preconditioner which depends on a
% parameter @var{k}.
%
% @example
% @group
% function y = apply_M (x, varargin)
% K = varargin@{1@};
% y = x;
% y(1:K) = x(1:K)./[1:K]';
%
%
% [x, flag, relres, iter, resvec, eigest] = ...
% pcg (A, b, [], [], 'apply_m', [], [], 3)
% @end group
% @end example
%
% @sc{References}
%
% [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations',
% SIAM, 1995 (the base PCG algorithm)
%
% [2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996
% (condition number estimate from PCG) Revised version of this book is
% available online at http://www-users.cs.umn.edu/~saad/books.html
%
%
% @seealso{sparse, pcr}
% @end deftypefn