— Function File: x = pcr (a, b, tol, maxit, m, x0, ...)
— Function File: [x, flag, relres, iter, resvec] = pcr (...)

Solves the linear system of equations a * x = b by means of the Preconditioned Conjugate Residuals iterative method. The input arguments are

• 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 a * x. In principle a should be symmetric and non-singular; if pcr finds a to be numerically singular, you will get a warning message and the flag output parameter will be set.
• b is the right hand side vector.
• tol is the required relative tolerance for the residual error, b - a * x. The iteration stops if norm (b - a * x) <= tol * norm (b - a * x0). If tol is empty or is omitted, the function sets tol = 1e-6 by default.
• maxit is the maximum allowable number of iterations; if [] is supplied for maxit, or pcr has less arguments, a default value equal to 20 is used.
• m is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by pcr P * x = m \ b, with P = m \ a. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrix m, the user may pass a function which returns the results of applying the inverse of m to a vector (usually this is the preferred way of using the preconditioner). If [] is supplied for m, or m is omitted, no preconditioning is applied.
• x0 is the initial guess. If x0 is empty or omitted, the function sets x0 to a zero vector by default.

The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (a or m) which are passed to pcr. See the examples below for further details. The output arguments are

• x is the computed approximation to the solution of a * x = b.
• flag reports on the convergence. flag = 0 means the solution converged and the tolerance criterion given by tol is satisfied. flag = 1 means that the maxit limit for the iteration count was reached. flag = 3 reports t pcr breakdown, see [1] for details.
• relres is the ratio of the final residual to its initial value, measured in the Euclidean norm.
• iter is the actual number of iterations performed.
• resvec describes the convergence history of the method, so that resvec (i) contains the Euclidean norms of the residualafter the (i-1)-th iteration, i = 1,2, ..., iter+1.

Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A)

          	N = 10;
          	A = diag([1:N]); A = sparse(A);
          	b = rand(N,1);
     

Example 1: Simplest use of pcr

            x = pcr(A, b)
     

Example 2: pcr with a function which computes a * x.

            function y = applyA(x)
              y = [1:10]'.*x;
            endfunction
          
            x = pcr('applyA',b)
     

Example 3: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix a is trivial) is defined as a function

            function y = applyM(x)
              K = floor(length(x)-2);
              y = x;
              y(1:K) = x(1:K)./[1:K]';
            endfunction
          
            [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM')
            semilogy([1:iter+1], resvec);
     

Example 4: Finally, a preconditioner which depends on a parameter k.

            function y = applyM(x, varargin)
              K = varargin{1};
              y = x; y(1:K) = x(1:K)./[1:K]';
            endfunction
          
            [x, flag, relres, iter, resvec] = pcr(A,b,[],[],'applyM',[],3)
     

References

[1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of Equations", section 9.5.4; Springer, 1994

See also: sparse, pcg.