Syntax

Description

x = gmres(A,b) attempts to solve the system of linear equations A*x = b for x. The n-by-n 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(b-A*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, 1e-6.

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 all-zero vector.

[x,flag] = gmres(A,b,...) also returns a convergence flag:

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(b-A*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.

Examples

Using gmres with a Matrix Input

A = gallery('wilk',21);  
b = sum(A,2);
tol = 1e-12;  
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.3e-013

Using gmres with a Function Handle

This example replaces the matrix A in the previous example with a handle to a matrix-vector 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 = 1e-12;  maxit = 15; 
x1 = gmres(@afun,b,10,tol,maxit,@mfun);
 
    function y = afun(x)
        y = [0; x(1:n-1)] + ...
              [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...
              [x(2:n); 0];
    end
 
    function y = mfun(r)
        y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/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.1e-013.

Using a Preconditioner without Restart

This example demonstrates the use of a preconditioner without restarting gmres.

Load west0479, a real 479-by-479 nonsymmetric sparse matrix.

load west0479;
A = west0479;

Set the tolerance and maximum number of iterations.

tol = 1e-12;
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 1e-12 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',1e-5));

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',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);

fl1 is 0 because gmres drives the relative residual to 9.5436e-14 (the value of rr1). The relative residual is less than the prescribed tolerance of 1e-12 at the sixth iteration (the value of it1(2)) when preconditioned by the incomplete LU factorization with a drop tolerance of 1e-6. The output, rv1(1), is norm(M\b), where M = L*U. The output, rv1(7), is norm(U\(L\(b-A*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');

Using a Preconditioner with Restart

This example demonstrates the use of a preconditioner with restarted gmres.

Load west0479, a real 479-by-479 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',1e-6));

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 = 1e-12;
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 1e-12.

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');

References

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. 856-869.

See Also

bicg | bicgstab | cgs | ilu | lsqr | minres | mldivide | pcg | qmr | symmlq