## Copyright (C) 2004-2018 Piotr Krzyzanowski
## Copyright (C) 2016-2018 Cristiano Dorigo
##
## 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
## .
## -*- texinfo -*-
## @deftypefn {} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{m1}, @var{m2}, @var{x0}, @dots{})
## @deftypefnx {} {@var{x} =} pcg (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, [], @var{x0}, @dots{})
## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}, @var{eigest}] =} pcg (@var{A}, @var{b}, @dots{})
##
## Solve the linear system of equations @w{@code{@var{A} * @var{x} = @var{b}}}
## by means of the Preconditioned Conjugate Gradient iterative method.
##
## The input arguments are:
##
## @itemize
## @item @var{A} is the matrix of the linear system and it must be square.
## @var{A} can be passed as a matrix, function handle, or inline function
## @code{Afun} such that @code{Afun(x) = A * x}. Additional parameters to
## @code{Afun} may be passed after @var{x0}.
##
## @var{A} has to be Hermitian and Positive Definite (@nospell{HPD})@. If
## @code{pcg} detects @var{A} not to be positive definite, a warning is printed
## and the @var{flag} output is set.
##
## @item
## @var{b} is the right-hand side vector.
##
## @item
## @var{tol} is the required relative tolerance for the residual error,
## @w{@code{@var{b} - @var{A} * @var{x}}}. The iteration stops if
## @w{@code{norm (@var{b} - @var{A} * @var{x})} @leq{}
## @w{@code{@var{tol} * norm (@var{b})}}}.
## If @var{tol} is omitted or empty, then a tolerance of 1e-6 is used.
##
## @item
## @var{maxit} is the maximum allowed number of iterations; if @var{maxit}
## is omitted or empty then a value of 20 is used.
##
## @item
## @var{m} is a @nospell{HPD} preconditioning matrix. For any decomposition
## @code{@var{m} = @var{p1} * @var{p2}} such that
## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2})}} is @nospell{HPD}, the
## conjugate gradient method is formally applied to the linear system
## @w{@code{inv (@var{p1}) * @var{A} * inv (@var{p2}) * @var{y} = inv
## (@var{p1}) * @var{b}}},
## with @code{@var{x} = inv (@var{p2}) * @var{y}} (split preconditioning).
## In practice, at each iteration of the conjugate gradient method a
## linear system with matrix @var{m} is solved with @code{mldivide}.
## If a particular factorization
## @code{@var{m} = @var{m1} * @var{m2}} is available (for instance, an
## incomplete Cholesky factorization of @var{a}), the two matrices
## @var{m1} and @var{m2} can be passed and the relative linear systems
## are solved with the @code{mldivide} operator.
## 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.
## If @var{m1} is omitted or empty @code{[]}, then no preconditioning
## is applied. If no factorization of @var{m} is available, @var{m2}
## can be omitted or left [], and the input variable @var{m1} can be
## used to pass the preconditioner @var{m}.
##
## @item
## @var{x0} is the initial guess. If @var{x0} is omitted or empty then 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
## an appropriate manner to any of the functions (@var{A} or @var{m1} or
## @var{m2}) that have been given 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
## @w{@code{@var{A} * @var{x} = @var{b}}}. If the algorithm did not converge,
## then @var{x} is the iteration which has the minimum residual.
##
## @item
## @var{flag} reports on the convergence:
##
## @itemize
## @item 0: The algorithm converged to within the prescribed tolerance.
##
## @item 1: The algorithm did not converge and it reached the maximum
## number of iterations.
##
## @item 2: The preconditioner matrix is singular.
##
## @item 3: The algorithm stagnated, i.e., the absolute value of the
## difference between the current iteration @var{x} and the previous is less
## than @code{@var{eps} * norm (@var{x},2)}.
##
## @item 4: The algorithm detects that the input (preconditioned) matrix is not
## @nospell{HPD}.
## @end itemize
##
## @item
## @var{relres} is the ratio of the final residual to its initial value,
## measured in the Euclidean norm.
##
## @item
## @var{iter} indicates the iteration of @var{x} which it was
## computed. Since the output @var{x} corresponds to the minimal
## residual solution, the total number of iterations that
## the method performed is given by @code{length(resvec) - 1}.
##
## @item
## @var{resvec} describes the convergence history of the method.
## @code{@var{resvec} (@var{i}, 1)} is the Euclidean norm of the residual, and
## @code{@var{resvec} (@var{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{@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
## @w{@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.
## @end itemize
##
##
## Let us consider a trivial problem with a tridiagonal matrix
##
## @example
## @group
## n = 10;
## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
## b = A * ones (n, 1);
## M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
## M2 = M1';
## M = M1 * M2;
## Afun = @@(x) A * x;
## Mfun = @@(x) M \ x;
## M1fun = @@(x) M1 \ x;
## M2fun = @@(x) M2 \ x;
## @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
## x = pcg (Afun, b)
## @end example
##
## @sc{Example 3:} @code{pcg} with a preconditioner matrix @var{M}
##
## @example
## x = pcg (A, b, 1e-06, 100, M)
## @end example
##
## @sc{Example 4:} @code{pcg} with a function as preconditioner
##
## @example
## x = pcg (Afun, b, 1e-6, 100, Mfun)
## @end example
##
## @sc{Example 5:} @code{pcg} with preconditioner matrices @var{M1}
## and @var{M2}
##
## @example
## x = pcg (A, b, 1e-6, 100, M1, M2)
## @end example
##
## @sc{Example 6:} @code{pcg} with functions as preconditioners
##
## @example
## x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)
## @end example
##
## @sc{Example 7:} @code{pcg} with as input a function requiring an argument
##
## @example
## @group
## function y = Ap (A, x, p) # compute A^p * x
## y = x;
## for i = 1:p
## y = A * y;
## endfor
## endfunction
## Apfun = @@(x, p) Ap (A, x, p);
## x = pcg (Apfun, b, [], [], [], [], [], 2);
## @end group
## @end example
##
## @sc{Example 8:} explicit example to show that @code{pcg} uses a
## split preconditioner
##
## @example
## @group
## M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
## M2 = M1';
## M = M1 * M2;
##
## ## reference solution computed by pcg after two iterations
## [x_ref, fl] = pcg (A, b, [], 2, M)
##
## ## split preconditioning
## [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
## x = M2 \ y # compare x and x_ref
##
## @end group
## @end example
##
## References:
##
## @enumerate
## @item
## C.T. Kelley, @cite{Iterative Methods for Linear and Nonlinear Equations},
## SIAM, 1995. (the base PCG algorithm)
##
## @item
## @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear Systems},
## @nospell{PWS} 1996. (condition number estimate from PCG)
## Revised version of this book is available online at
## @url{https://www-users.cs.umn.edu/~saad/books.html}
## @end enumerate
##
## @seealso{sparse, pcr, gmres, bicg, bicgstab, cgs}
## @end deftypefn
## Author: Piotr Krzyzanowski
## Modified by: Vittoria Rezzonico
## - Add the ability to provide the pre-conditioner as two separate matrices
function [x_min, flag, relres, iter_min, resvec, eigest] =...
bpcg (A, b, tol = [], maxit = [], M1 = [], M2 = [], x0 = [], varargin)
## Insert the default input (if necessary)
[tol, maxit, x0] = __default__input__ ({1e-6, min(rows (b), 20),...
zeros(size (b))}, tol, maxit, x0);
if (tol >= 1)
warning ("Input tol is bigger than 1. \n Try to use a smaller tolerance.");
elseif (tol <= eps / 2)
warning ("Input tol may not be achievable by pcg. \n Try to use a bigger tolerance");
endif
## Check if the input data A,b,m1,m2 are consistent (i.e. if they are
## matrix or function handle)
[Afun, M1fun, M2fun] = __alltohandles__ (A, b, M1, M2, "pcg");
maxit += 2;
n_arg_out = nargout;
## Set Initial data
b_norm = norm (b,'fro');
if (b_norm == 0)
if (n_arg_out < 2)
printf("The right hand side vector is all zero so pcg \n");
printf ("returned an all zero solution without iterating.\n");
endif
x_min = b;
flag = 0;
relres = 0;
resvec = 0;
iter_min = 0;
eigest = [NaN, NaN];
return
endif
x = x_pr = x_min = x0;
## x_pr (x previous) needs to check the stagnation
## x_min needs to save the iterated with minimum residual
r = b - feval (Afun, x, varargin{:});
iter = 2;
iter_min = 0;
flag = 1;
resvec = zeros (maxit + 1, 2);
resvec(1, 1) = norm (r,'fro');
p = zeros (size (b));
alpha = den = tau = 1;
if (n_arg_out > 5)
ncase = size (b, 2);
T = zeros (maxit*ncase, maxit*ncase);
else
T = [];
endif
while (resvec(iter-1,1) > tol * b_norm && iter < maxit)
if (iter == 2) # Check whether M1 or M2 are singular
try
warning ("error","Octave:singular-matrix","local")
z = feval (M1fun, r, varargin{:});
z = feval (M2fun, z, varargin{:});
catch
flag = 2;
break;
end_try_catch
else
z = feval (M1fun, r, varargin{:});
z = feval (M2fun, z, varargin{:});
endif
old_tau = tau;
tau = z' * r;
resvec(iter - 1, 2) = norm (tau^(1/2), 'fro');
beta = pinv(old_tau) * tau;
p = z + p * beta;
w = feval (Afun, p, varargin{:});
## Needed only for eigest.
old_den = den;
old_alpha = alpha;
den = p' * w;
## Diagonalize the small matrices den and tau and pinvert den
ttl = eig(tau);
tdl = eig(den);
alpha = pinv(den)*tau;
## Check if alpha is negative and/or if it has a consistent
## imaginary part: if yes then A probably is not positive definite
if ((max (abs (imag (ttl))) >= max (abs (real (ttl))) * tol) || ...
(min (real (ttl)) <= 0) || ...
(max (abs (imag (tdl))) >= max (abs (real (tdl))) * tol) || ...
(min (real (tdl)) <= 0))
flag = 4;
break;
endif
x += p * alpha;
r -= w * alpha;
resvec(iter, 1) = norm (r, 'fro');
## Chek if the iterated has minimum residual
if (resvec (iter,1) <= resvec (iter_min + 1,1))
x_min = x;
iter_min = iter - 1;
endif
if (n_arg_out > 5 && iter > 2)
# Build the block-tridiagonal matrix
T11 = pinv(old_tau^(1/2))*old_den*pinv(old_tau^(1/2));
T12 = pinv(old_tau^(1/2))*old_den*pinv(old_tau)*tau^(1/2);
T22 = pinv(tau^(1/2))*(beta'*old_den*beta)*pinv(tau^(1/2));
midm1 = ncase*(iter-2)+1;
mid = ncase*iter;
T(midm1:mid, midm1:mid) = T(midm1:mid, midm1:mid) + [T11,T12;T12',T22];
endif
iter += 1;
if (norm (x - x_pr) <= eps * norm (x)) # Check the stagnation
flag = 3;
break;
endif
x_pr = x;
endwhile
if (n_arg_out > 5)
## Apply the preconditioner once more and finish with the precond
## residual.
z = feval (M1fun, r, varargin{:});
z = feval (M2fun, z, varargin{:});
endif
## (Eventually) computes the eigenvalue of inv(m2)*inv(m1)*A
if (n_arg_out > 5)
if (flag != 4)
if (iter > 3)
T = T(1+ncase:end-2*(ncase), 1+ncase:end-2*(ncase) );
l = eig (T);
eigest = [min(l), max(l)];
else
eigest = [NaN, NaN];
warning ("pcg: eigenvalue estimate failed: iteration converged too fast");
endif
else
eigest = [NaN, NaN];
warning ('pcg: eigenvalue estimate failed: matrix not positive definite?')
endif
resvec(iter - 1, 2) = norm ((r' * z)^(1/2), 'fro');
resvec = resvec (1:(iter-1), :);
else
eigest = [NaN, NaN];
resvec = resvec(1:(iter-1),1);
endif
## Set the last variables
if (flag == 2)
relres = 1;
elseif (resvec (1, 1) == 0)
relres = 0;
else
relres = resvec(iter_min+1, 1) ./ resvec(1, 1);
endif
iter -= 2; # compatibility
## Set the flag in the proper way if flag not 3, 4 or 2
if (flag == 2)
flag = 2;
elseif (flag == 1) && (relres <= tol)
flag = 0;
endif
if (n_arg_out < 2)
switch (flag)
case {0}
printf ("pcg converged at iteration %d ", iter_min);
printf ("with relative residual %d\n", relres);
case {1}
printf ("pcg stopped at iteration %d ", iter+1);
printf ("without converging to the desired tolerance %d ", tol);
printf ("because the maximum number of iteration was reached, \n");
printf ("The iterated returned (number %d) ",iter_min);
printf ("has relative residual %d \n", relres);
case {2}
printf ("pcg stopped at iteration %d ", iter+1)
printf ("without converging to the desired tolerance %d ", tol);
printf ("because the preconditioned matrix is singular.\n");
printf ("The iterated returned (number %d) ", iter_min);
printf ("has relative residual %d \n", relres);
case {3}
printf ("pcg stopped at iteration %d ", iter+1);
printf ("without converging to the desired tolerance %d ", tol);
printf ("because of stagnation. \n");
printf ("The iterated returned (number %d) ", iter_min);
printf ("has relative residual %d.\n", relres);
case {4}
printf ("pcg stopped at iteration %d ", iter + 1);
printf ("without converging to the desired tolerance %d ",tol);
printf ("because the (preconditioned) matrix is not positive definite. \n");
printf ("The iterate returned (number %d) ", iter_min);
printf ("has relative residual %d \n", relres);
endswitch
endif
endfunction
%!demo # simplest use
%! n = 10;
%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
%! b = A * ones (n, 1);
%! M1 = ichol (A); # for this tridiagonal case it corresponds to chol (A)'
%! M2 = M1';
%! M = M1 * M2;
%! x = pcg (A, b);
%! Afun = @(x) A * x;
%! x = pcg (Afun, b);
%! x = pcg (A, b, 1e-6, 100, M);
%! x = pcg (A, b, 1e-6, 100, M1, M2);
%! Mfun = @(x) M \ x;
%! x = pcg (Afun, b, 1e-6, 100, Mfun);
%! M1fun = @(x) M1 \ x;
%! M2fun = @(x) M2 \ x;
%! x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun);
%! function y = Ap (A, x, p) # compute A^p * x
%! y = x;
%! for i = 1:p
%! y = A * y;
%! endfor
%! endfunction
%! Afun = @(x, p) Ap (A, x, p);
%! ## solution of A^2 * x = b
%! x = pcg (Afun, b, [], [], [], [], [], 2);
%!demo
%! n = 10;
%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
%! b = A * ones (n, 1);
%! M1 = ichol (A + 0.1 * eye (n)); # Perturb the factorization of A
%! M2 = M1';
%! M = M1 * M2;
%!
%! ## Reference solution computed by pcg after two iterations
%! [x_ref, fl] = pcg (A, b, [], 2, M);
%! x_ref
%!
%! ## Split preconditioning
%! [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2);
%! x = M2 \ y # compare x and x_ref
%!test
%! ## Check that all type of inputs work
%! A = toeplitz (sparse ([2, 1 ,0, 0, 0]));
%! b = A * ones (5, 1);
%! M1 = diag (sqrt (diag (A)));
%! M2 = M1; # M1 * M2 is the Jacobi preconditioner
%! Afun = @(z) A*z;
%! M1_fun = @(z) M1 \ z;
%! M2_fun = @(z) M2 \ z;
%! [x, flag, ~, iter] = pcg (A,b);
%! assert (flag, 0);
%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1 * M2);
%! assert (flag, 0);
%! [x, flag, ~ , iter] = pcg (A, b, [], [], M1, M2);
%! assert (flag, 0);
%! [x, flag] = pcg (A, b, [], [], M1_fun, M2_fun);
%! assert (flag, 0);
%! [x, flag] = pcg (A, b,[],[], M1_fun, M2);
%! assert (flag, 0);
%! [x, flag] = pcg (A, b,[],[], M1, M2_fun);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b,[],[], M1 * M2);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b,[],[], M1, M2);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b,[],[], M1, M2_fun);
%! assert (flag, 0);
%! [x, flag] = pcg (Afun, b,[],[], M1_fun, M2_fun);
%! assert (flag, 0);
%!test
%! ## solve a small diagonal system
%! N = 10;
%! A = diag ([1:N]); b = rand (N, 1);
%! [x, flag] = pcg (A, b, [], N+1);
%! assert (flag, 0);
%! assert (norm (b - A*x) / norm (b), 0, 1e-6);
%!test
%! ## A is not positive definite
%! ## The indefiniteness of A is detected.
%! N = 10;
%! A = -diag ([1:N]); b = sum (A, 2);
%! [x, flag] = pcg (A, b, [], N + 1);
%! assert (flag, 4);
%!test
%! ## solve tridiagonal system, do not converge in default 20 iterations
%! N = 100;
%! A = zeros (N, N);
%! for i = 1 : N - 1 # form 1-D Laplacian matrix
%! A(i:i+1, i:i+1) = [2 -1; -1 2];
%! endfor
%! b = ones (N, 1);
%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, 1e-12);
%! assert (flag);
%! assert (relres >= 1.0);
%!warning
%! ## solve tridiagonal system with "perfect" preconditioner which converges
%! ## in one iteration, so the eigest does not work and issues a warning.
%! N = 100;
%! A = zeros (N, N);
%! for i = 1 : N - 1 # form 1-D Laplacian matrix
%! A(i:i+1, i:i+1) = [2 -1; -1 2];
%! endfor
%! b = ones (N, 1);
%! [x, flag, relres, iter, resvec, eigest] = pcg (A, b, [], [], A, [], b);
%! assert (flag, 0);
%! assert (norm (b - A*x) / norm (b), 0, 1e-6);
%!
%! assert (isnan (eigest), isnan ([NaN, NaN]));
%!test
%! ## pcg detect a non-Hermitian matrix, with a considerable imaginary part.
%! ## In this example, Matlab does not recognize the wrong type of matrix and
%! ## makes iterations until it reaches maxit.
%! N = 10;
%! A = diag (1:N) + 1e-4*i;
%! b = ones (N, 1);
%! [x, flag] = pcg (A, b, []);
%! assert (flag, 4);
%!test
%! ## The imaginary part is not influent (it is too small), so pcg doesn't stop
%! N = 10;
%! A = diag (1:N) + 1e-10*i;
%! b = ones (N, 1);
%! [x, flag] = pcg (A, b, [], N+1);
%! assert (flag, 0);
%! assert (norm (b - A*x) / norm (b), 0, 1e-6);
%!test
%! ## pcg solves linear system with A Hermitian positive definite
%! N = 20;
%! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
%! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)]));
%! b = A * ones (N, 1);
%! Hermitian_A = ishermitian (A);
%! [x, flag] = pcg (A, b, [], 2*N);
%! assert (Hermitian_A, true);
%! assert (flag, 0);
%! assert (x, ones (N, 1), -1e-4);
%!testif HAVE_CHOLMOD
%! ## pcg solves preconditioned linear system with A HPD
%! N = 20;
%! A = sparse (toeplitz ([4, 1, zeros(1, 18)])) + ...
%! i * sparse (toeplitz ([0, 1, zeros(1, 18)], [0, -1, zeros(1,18)]));
%! b = A * ones (N, 1);
%! M2 = chol (A + 0.1 * eye (N)); # Factor of a perturbed matrix
%! M = M2' * M2;
%! Hermitian_A = ishermitian (A);
%! Hermitian_M = ishermitian (M);
%! [x, flag] = pcg (A, b, [], 2*N, M);
%! assert (Hermitian_A, true);
%! assert (Hermitian_M, true);
%! assert (flag, 0);
%! assert (x, ones (N, 1), -1e-4);
%!test
%! ## pcg recognizes that the preconditioner matrix is singular
%! N = 3;
%! A = toeplitz ([2, 1, 0]);
%! M = [1 0 0; 0 1 0; 0 0 0]; # the last row is zero
%! [x, flag] = pcg (A, ones (3, 1), [], [], M);
%! assert (flag, 2);
%!test
%! A = rand (4);
%! A = A' * A;
%! [x, flag] = pcg (A, zeros (4, 1), [], [], [], [], ones (4, 1));
%! assert (x, zeros (4, 1));
## Test return types
%!test
%! A = single (1);
%! b = 1;
%! [x, flag] = pcg (A, b);
%! assert (class (x), "single");
%!test
%! A = 1;
%! b = single (1);
%! [x, flag] = pcg (A, b);
%! assert (class (x), "single");
%!test
%! A = single (1);
%! b = single (1);
%! [x, flag] = pcg (A, b);
%! assert (class (x), "single");
%!test
%!function y = Afun (x)
%! A = toeplitz ([2, 1, 0, 0]);
%! y = A * x;
%!endfunction
%! [x, flag] = pcg ("Afun", [3; 4; 4; 3]);
%! assert (x, ones (4, 1), 1e-6);
%!test
%! ## unpreconditioned residual
%! A = toeplitz (sparse ([4, 1, 0, 0, 0]));
%! b = sum (A, 2);
%! M = toeplitz (sparse ([2, 1, 0, 0, 0]));
%! [x, flag, relres] = pcg (A, b, [], 2, M);
%! assert (norm (b - A * x) / norm (b), relres, 8 * eps);