Generalized minimum residual method (with restarts)
x = gmres(A,b)
gmres(A,b,restart)
gmres(A,b,restart,tol)
gmres(A,b,restart,tol,maxit)
gmres(A,b,restart,tol,maxit,M)
gmres(A,b,restart,tol,maxit,M1,M2)
gmres(A,b,restart,tol,maxit,M1,M2,x0)
[x,flag] = gmres(A,b,...
)
[x,flag,relres] = gmres(A,b,...
)
[x,flag,relres,iter] = gmres(A,b,...
)
[x,flag,relres,iter,resvec] = gmres(A,b,...
)
x = gmres(A,b)
attempts
to solve the system of linear equations A*x = b
for x
.
The n
byn
coefficient matrix A
must
be square and should be large and sparse. The column vector b
must
have length n
. A
can be a function
handle, afun
, such that afun(x)
returns A*x
.
For this syntax, gmres
does not restart; the maximum
number of iterations is min(n,10)
.
Parameterizing Functions explains how to provide additional
parameters to the function afun
, as well as the
preconditioner function mfun
described below,
if necessary.
If gmres
converges, a message to that effect
is displayed. If gmres
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.
gmres(A,b,restart)
restarts
the method every restart
inner iterations. The
maximum number of outer iterations is min(n/restart,10)
.
The maximum number of total iterations is restart*min(n/restart,10)
.
If restart
is n
or []
,
then gmres
does not restart and the maximum number
of total iterations is min(n,10)
.
gmres(A,b,restart,tol)
specifies
the tolerance of the method. If tol
is []
,
then gmres
uses the default, 1e6
.
gmres(A,b,restart,tol,maxit)
specifies
the maximum number of outer iterations, i.e., the total number of
iterations does not exceed restart*maxit
. If maxit
is []
then gmres
uses
the default, min(n/restart,10)
. If restart
is n
or []
,
then the maximum number of total iterations is maxit
(instead
of restart*maxit
).
gmres(A,b,restart,tol,maxit,M)
and
gmres(A,b,restart,tol,maxit,M1,M2)
use preconditioner M
or M
= M1*M2
and effectively solve the system inv(M)*A*x
= inv(M)*b
for x
. If M
is []
then gmres
applies
no preconditioner. M
can be a function handle mfun
such
that mfun(x)
returns M\x
.
gmres(A,b,restart,tol,maxit,M1,M2,x0)
specifies
the first initial guess. If x0
is []
,
then gmres
uses the default, an allzero vector.
[x,flag] = gmres(A,b,
also
returns a convergence flag:...
)




 Preconditioner 


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] = gmres(A,b,
also
returns the relative residual ...
)norm(bA*x)/norm(b)
.
If flag
is 0
, relres
<= tol
. The third output, relres
,
is the relative residual of the preconditioned system.
[x,flag,relres,iter] = gmres(A,b,
also
returns both the outer and inner iteration numbers at which ...
)x
was
computed, where 0 <= iter(1) <= maxit
and 0
<= iter(2) <= restart
.
[x,flag,relres,iter,resvec] = gmres(A,b,
also
returns a vector of the residual norms at each inner iteration. These
are the residual norms for the preconditioned system....
)
A = gallery('wilk',21); b = sum(A,2); tol = 1e12; maxit = 15; M1 = diag([10:1:1 1 1:10]); x = gmres(A,b,10,tol,maxit,M1);
displays the following message:
gmres(10) converged at outer iteration 2 (inner iteration 9) to a solution with relative residual 3.3e013
This example replaces the matrix A
in the
previous example with a handle to a matrixvector product function afun
,
and the preconditioner M1
with a handle to a backsolve
function mfun
. The example is contained in a
function run_gmres
that
Calls gmres
with the function
handle @afun
as its first argument.
Contains afun
and mfun
as
nested functions, so that all variables in run_gmres
are
available to afun
and mfun
.
The following shows the code for run_gmres
:
function x1 = run_gmres n = 21; b = afun(ones(n,1)); tol = 1e12; maxit = 15; x1 = gmres(@afun,b,10,tol,maxit,@mfun); function y = afun(x) y = [0; x(1:n1)] + ... [((n1)/2:1:0)'; (1:(n1)/2)'].*x + ... [x(2:n); 0]; end function y = mfun(r) y = r ./ [((n1)/2:1:1)'; 1; (1:(n1)/2)']; end end
When you enter
x1 = run_gmres;
MATLAB^{®} software displays the message
gmres(10) converged at outer iteration 2 (inner iteration 10) to a solution with relative residual 1.1e013.
This example demonstrates the use of a preconditioner without restarting gmres
.
Load west0479
, a real 479by479 nonsymmetric sparse matrix.
load west0479;
A = west0479;
Set the tolerance and maximum number of iterations.
tol = 1e12; maxit = 20;
Define b
so that the true solution is a vector of all ones.
b = full(sum(A,2)); [x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0
is 1 because gmres
does not converge to the requested tolerance 1e12
within the requested 20 iterations. The best approximate solution that gmres
returns is the last one (as indicated by it0(2) = 20
). MATLAB stores the residual history in rv0
.
Plot the behavior of gmres
.
semilogy(0:maxit,rv0/norm(b),'o'); xlabel('Iteration number'); ylabel('Relative residual');
The plot shows that the solution converges slowly. A preconditioner may improve the outcome.
Use ilu
to form the preconditioner, since A
is nonsymmetric.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e5));
Error using ilu There is a pivot equal to zero. Consider decreasing the drop tolerance or consider using the 'udiag' option.
Note MATLAB cannot construct the incomplete LU as it would result in a singular factor, which is useless as a preconditioner.
As indicated by the error message, try again with a reduced drop tolerance.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e6)); [x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1
is 0 because gmres
drives the relative residual to 9.5436e14
(the value of rr1
). The relative residual is less than the prescribed tolerance of 1e12
at the sixth iteration (the value of it1(2)
) when preconditioned by the incomplete LU factorization with a drop tolerance of 1e6
. The output, rv1(1)
, is norm(M\b)
, where M = L*U
. The output, rv1(7)
, is norm(U\(L\(bA*x1)))
.
Follow the progress of gmres
by plotting the relative residuals at each iteration starting from the initial estimate (iterate number 0).
semilogy(0:it1(2),rv1/norm(b),'o'); xlabel('Iteration number'); ylabel('Relative residual');
This example demonstrates the use of a preconditioner with restarted gmres
.
Load west0479
, a real 479by479 nonsymmetric sparse matrix.
load west0479;
A = west0479;
Define b
so that the true solution is a vector of all ones.
b = full(sum(A,2));
Construct an incomplete LU preconditioner as in the previous example.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e6));
The benefit to using restarted gmres
is to limit the amount of memory required to execute the method. Without restart, gmres
requires maxit
vectors of storage to keep the basis of the Krylov subspace. Also, gmres
must orthogonalize against all of the previous vectors at each step. Restarting limits the amount of workspace used and the amount of work done per outer iteration. Note that even though preconditioned gmres
converged in six iterations above, the algorithm allowed for as many as twenty basis vectors and therefore, allocated all of that space up front.
Execute gmres(3)
, gmres(4)
, and gmres(5)
tol = 1e12; maxit = 20; re3 = 3; [x3,fl3,rr3,it3,rv3] = gmres(A,b,re3,tol,maxit,L,U); re4 = 4; [x4,fl4,rr4,it4,rv4] = gmres(A,b,re4,tol,maxit,L,U); re5 = 5; [x5,fl5,rr5,it5,rv5] = gmres(A,b,re5,tol,maxit,L,U);
fl3
, fl4
, and fl5
are all 0 because in each case restarted gmres
drives the relative residual to less than the prescribed tolerance of 1e12
.
The following plots show the convergence histories of each restarted gmres
method. gmres(3)
converges at outer iteration 5, inner iteration 3 (it3 = [5, 3]
) which would be the same as outer iteration 6, inner iteration 0, hence the marking of 6 on the final tick mark.
figure semilogy(1:1/3:6,rv3/norm(b),'o'); h1 = gca; h1.XTick = [1:1/3:6]; h1.XTickLabel = ['1';' ';' ';'2';' ';' ';'3';' ';' ';'4';' ';' ';'5';' ';' ';'6';]; title('gmres(3)') xlabel('Iteration number'); ylabel('Relative residual'); figure semilogy(1:1/4:3,rv4/norm(b),'o'); h2 = gca; h2.XTick = [1:1/4:3]; h2.XTickLabel = ['1';' ';' ';' ';'2';' ';' ';' ';'3']; title('gmres(4)') xlabel('Iteration number'); ylabel('Relative residual'); figure semilogy(1:1/5:2.8,rv5/norm(b),'o'); h3 = gca; h3.XTick = [1:1/5:2.8]; h3.XTickLabel = ['1';' ';' ';' ';' ';'2';' ';' ';' ';' ']; title('gmres(5)') xlabel('Iteration number'); ylabel('Relative residual');
Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.
Saad, Youcef and Martin H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., July 1986, Vol. 7, No. 3, pp. 856869.