Preconditioned conjugate gradients method
x = pcg(A,b)
pcg(A,b,tol)
pcg(A,b,tol,maxit)
pcg(A,b,tol,maxit,M)
pcg(A,b,tol,maxit,M1,M2)
pcg(A,b,tol,maxit,M1,M2,x0)
[x,flag] = pcg(A,b,...)
[x,flag,relres] = pcg(A,b,...)
[x,flag,relres,iter] = pcg(A,b,...)
[x,flag,relres,iter,resvec] = pcg(A,b,...)
x = pcg(A,b)
attempts to
solve the system of linear equations A*x=b
for x
.
The n
byn
coefficient matrix A
must
be symmetric and positive definite, and should also be large and sparse.
The column vector b
must have length n
.
You also can specify A
to be a function handle, afun
,
such that afun(x)
returns A*x
.
Parameterizing Functions explains how to provide additional
parameters to the function afun
, as well as the
preconditioner function mfun
described below, if
necessary.
If pcg
converges, a message to that effect
is displayed. If pcg
fails to converge after the
maximum number of iterations or halts for any reason, a warning message
is printed displaying the relative residual norm(bA*x)/norm(b)
and
the iteration number at which the method stopped or failed.
pcg(A,b,tol)
specifies
the tolerance of the method. If tol
is []
,
then pcg
uses the default, 1e6
.
pcg(A,b,tol,maxit)
specifies
the maximum number of iterations. If maxit
is []
,
then pcg
uses the default, min(n,20)
.
pcg(A,b,tol,maxit,M)
and pcg(A,b,tol,maxit,M1,M2)
use symmetric positive definite preconditioner M
or M
= M1*M2
and effectively solve the system inv(M)*A*x
= inv(M)*b
for x
. If M
is []
then pcg
applies
no preconditioner. M
can be a function handle mfun
such
that mfun(x)
returns M\x
.
pcg(A,b,tol,maxit,M1,M2,x0)
specifies
the initial guess. If x0
is []
,
then pcg
uses the default, an allzero vector.
[x,flag] = pcg(A,b,...)
also
returns a convergence flag.
Flag  Convergence 





 Preconditioner 


 One of the scalar quantities calculated during 
Whenever flag
is not 0
,
the solution x
returned is that with minimal norm
residual computed over all the iterations. No messages are displayed
if the flag
output is specified.
[x,flag,relres] = pcg(A,b,...)
also
returns the relative residual norm(bA*x)/norm(b)
.
If flag
is 0
, relres
<= tol
.
[x,flag,relres,iter] = pcg(A,b,...)
also
returns the iteration number at which x
was computed,
where 0 <= iter <= maxit
.
[x,flag,relres,iter,resvec] = pcg(A,b,...)
also
returns a vector of the residual norms at each iteration including norm(bA*x0)
.
This example shows how to use pcg with a matrix input and with a function handle.
n1 = 21; A = gallery('moler',n1); b1 = sum(A,2); tol = 1e6; maxit = 15; M1 = spdiags((1:n1)',0,n1,n1); [x1,flag1,rr1,iter1,rv1] = pcg(A,b1,tol,maxit,M1);
Alternatively, you can use the following function in place of
the matrix A
:
function y = applyMoler(x) y = x; y(end1:1:1) = y(end1:1:1)  cumsum(y(end:1:2)); y(2:end) = y(2:end)  cumsum(y(1:end1));
By using this function, you can solve larger systems more efficiently
as there is no need to store the entire matrix A
:
n2 = 21; b2 = applyMoler(ones(n2,1)); tol = 1e6; maxit = 15; M2 = spdiags((1:n2)',0,n2,n2); [x2,flag2,rr2,iter2,rv2] = pcg(@applyMoler,b2,tol,maxit,M2);
pcg
with a PreconditionerThis example demonstrates how to use a preconditioner matrix with pcg
.
Create an input matrix and try to solve the system with pcg
.
A = delsq(numgrid('S',100));
b = ones(size(A,1),1);
[x0,fl0,rr0,it0,rv0] = pcg(A,b,1e8,100);
fl0
is 1
because pcg
does not converge to the requested tolerance of 1e8
within the requested maximum 100 iterations. A preconditioner can make the system converge more quickly.
Use ichol
with only one input argument to construct an incomplete Cholesky factorization with zero fill.
L = ichol(A); [x1,fl1,rr1,it1,rv1] = pcg(A,b,1e8,100,L,L');
fl1
is 0
because pcg
drives the relative residual to 9.8e09
(the value of rr1
) which is less than the requested tolerance of 1e8
at the seventyseventh iteration (the value of it1
) when preconditioned by the zerofill incomplete Cholesky factorization. rv1(1) = norm(b)
and rv1(78) = norm(bA*x1)
.
The previous matrix represents the discretization of the Laplacian on a 100x100 grid with Dirichlet boundary conditions. This means that a modified incomplete Cholesky preconditioner might perform even better.
Use the michol
option to create a modified incomplete Cholesky preconditioner.
L = ichol(A,struct('michol','on')); [x2,fl2,rr2,it2,rv2] = pcg(A,b,1e8,100,L,L');
In this case you attain convergence in only fortyseven iterations.
You can see how the preconditioners affect the rate of convergence of pcg
by plotting each of the residual histories starting from the initial estimate (iterate number 0
).
figure; semilogy(0:it0,rv0/norm(b),'b.'); hold on; semilogy(0:it1,rv1/norm(b),'r.'); semilogy(0:it2,rv2/norm(b),'k.'); legend('No Preconditioner','IC(0)','MIC(0)'); xlabel('iteration number'); ylabel('relative residual'); hold off;
[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.