Function Reference
— 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 residual after 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

Demonstration 1

The following code


	# Simplest usage of PCR (see also 'help pcr')

	N = 20; 
	A = diag(linspace(-3.1,3,N)); b = rand(N,1); y = A\b; #y is the true solution
  	x = pcr(A,b);
	printf('The solution relative error is %g\n', norm(x-y)/norm(y));

	# You shouldn't be afraid if PCR issues some warning messages in this
	# example: watch out in the second example, why it takes N iterations 
	# of PCR to converge to (a very accurate, by the way) solution

Produces the following output

warning: pcr: maximum number of iterations (20) reached

warning: the initial residual norm was reduced 2.41124e+10 times.

The solution relative error is 8.68184e-12

Demonstration 2

The following code


	# Full output from PCR
	# We use this output to plot the convergence history  

	N = 20; 
	A = diag(linspace(-3.1,30,N)); b = rand(N,1); X = A\b; #X is the true solution
  	[x, flag, relres, iter, resvec] = pcr(A,b);
	printf('The solution relative error is %g\n', norm(x-X)/norm(X));
	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||/||b||)');
	semilogy([0:iter],resvec/resvec(1),'o-g;relative residual;');

Produces the following output

The solution relative error is 9.86986e-14

and the following figure

Demonstration 3

The following code


	# Full output from PCR
	# We use indefinite matrix based on the Hilbert matrix, with one
	# strongly negative eigenvalue
	# Hilbert matrix is extremely ill conditioned, so is ours, 
	# and that's why PCR WILL have problems

	N = 10; 
	A = hilb(N); A(1,1)=-A(1,1); b = rand(N,1); X = A\b; #X is the true solution
	printf('Condition number of A is   %g\n', cond(A));
  	[x, flag, relres, iter, resvec] = pcr(A,b,[],200);
	if (flag == 3)
	  printf('PCR breakdown. System matrix is [close to] singular\n');
	end
	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
	semilogy([0:iter],resvec,'o-g;absolute residual;');

Produces the following output

Condition number of A is   8.64683e+12
PCR breakdown. System matrix is [close to] singular

and the following figure

Demonstration 4

The following code


	# Full output from PCR
	# We use an indefinite matrix based on the 1-D Laplacian matrix for A, 
	# and here we have cond(A) = O(N^2)
	# That's the reason we need some preconditioner; here we take
	# a very simple and not powerful Jacobi preconditioner, 
	# which is the diagonal of A

	# Note that we use here indefinite preconditioners!

	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
	A = [A, zeros(size(A)); zeros(size(A)), -A];
	b = rand(2*N,1); X = A\b; #X is the true solution
	maxit = 80;
	printf('System condition number is %g\n',cond(A));
	# No preconditioner: the convergence is very slow!

  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit);
	title('Convergence history'); xlabel('Iteration'); ylabel('log(||b-Ax||)');
	semilogy([0:iter],resvec,'o-g;NO preconditioning: absolute residual;');

	pause(1);
	# Test Jacobi preconditioner: it will not help much!!!

	M = diag(diag(A)); # Jacobi preconditioner
  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
	hold on;
	semilogy([0:iter],resvec,'o-r;JACOBI preconditioner: absolute residual;');

	pause(1);
	# Test nonoverlapping block Jacobi preconditioner: this one should give
	# some convergence speedup!

	M = zeros(N,N);k=4;
	for i=1:k:N # get k x k diagonal blocks of A
		M(i:i+k-1,i:i+k-1) = A(i:i+k-1,i:i+k-1);
	endfor
	M = [M, zeros(size(M)); zeros(size(M)), -M];
  	[x, flag, relres, iter, resvec] = pcr(A,b,[],maxit,M);
	semilogy([0:iter],resvec,'o-b;BLOCK JACOBI preconditioner: absolute residual;');
	hold off;

Produces the following output

System condition number is 4133.64

and the following figure