MATLAB: code for restarted gmres

Question

I have a question about Matlab and restarted gmres. I would like to use gmres.m provided here. This code seems to be popular for the scientific computation newcomer. I also want to modify this code for my particular case.

But I have a question about recording the number of matrix-vector products and the 2-norm of each residual vector. If you use the built-in function gmres.m (e.g., gmres(30)) to solve a linear system, one obtains

iter = [17,20], % the total number of matrix-vector products is (17 - 1)*30 + 12 (= 492).

and resvec, a $493\times1$ vector.

Moreover, in the gmres.m function these two outputs are not available. As you know, we always call the gmres.m in MATLAB in the following style:

[x, flag, relres, iter, resvec] = gmres(A, b, restart, tol, maxit);

Then, I want to use these outputs to plot the convergence histories, i.e.,

X-axis: the number of matrix-vector products,

Y-axis: 2-norm of relative residual, i.e., norm(b - Ax_m)/norm(b);

So, how can I record this two outputs reasonably in the modified gmres.m? Despite my many attempts to deal with this problem, I still fail to manage it. It makes me very confused, especially for recording resvec. Can anyone provide suggestions for handling this problem?

Answer 1

Technically, the template Matlab code provided at Netlib already calculates the obtained residual at every iteration. It's just not recorded in the way that this is an output of the subroutine.

After each iteration iter the residual r is calculated, so one just needs to record it. I will do it right after it is calculated in the residual vector resvec(iter).

Now, the provided GMRES code is not optimal (for example, the residual r is calculated both before and after the loop over i-variable: construction of the orthonormal basis), so I am not going to optimize that behaviour, but just count matrix-vector products (MVP) every time they happen in mvp_count. Then, at the end of each iteration, I will record the number of actually performed MVPs in mvpvec(iter).

Both resvec and mvpvec are added to the output parameters of the gmres function which should provide enough information to do the plot.

Notes:

• mvpvec(iter) counts MVPs performed BY the end of the iterth iteration by this particular implementation of gmres.
• no normalization by the 2-norm of the RHS b is done (to avoid confusion between residul, relative error, etc.). Should be easily changed either inside or outside of the function.
• only MVPs with the supposedly large original matrix $A$ are being counted
• the initial residual and the corresponding MVP (before iterative process starts) are not recorded.

Here is the slightly modified Matlab code.

function [x, error, iter, flag, resvec, mvpvec] = gmres( A, x, b, M, restrt, max_it, tol )

%  -- Iterative template routine --
%     Univ. of Tennessee and Oak Ridge National Laboratory
%     October 1, 1993
%     Details of this algorithm are described in "Templates for the
%     Solution of Linear Systems: Building Blocks for Iterative
%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
%     Sligtly modified by Anton Menshov to answer CompSci SE question in 2018.
%
% [x, error, iter, flag, resvec, mvpvec] = gmres( A, x, b, M, restrt, max_it, tol )
%
% gmres.m solves the linear system Ax=b
% using the Generalized Minimal residual ( GMRESm ) method with restarts .
%
% input   A        REAL nonsymmetric positive definite matrix
%         x        REAL initial guess vector
%         b        REAL right hand side vector
%         M        REAL preconditioner matrix
%         restrt   INTEGER number of iterations between restarts
%         max_it   INTEGER maximum number of iterations
%         tol      REAL error tolerance
%
% output  x        REAL solution vector
%         error    REAL error norm
%         iter     INTEGER number of iterations performed
%         flag     INTEGER: 0 = solution found to tolerance
%                           1 = no convergence given max_it
%         resvec   REAL (iter x 1) vector of norm2 solution resuduals at each iteration
%         mvpvec   INTEGER (iter x 1) vector of ACTUALLY computed MVPs at each iteration

iter = 0;                                         % initialization
flag = 0;
mvp_count =0;

bnrm2 = norm( b );
if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

r = M \ ( b-A*x );     % not recorded and MVP is not counted
error = norm( r ) / bnrm2;
if ( error < tol ) return, end

[n,n] = size(A);                                  % initialize workspace
m = restrt;
V(1:n,1:m+1) = zeros(n,m+1);
H(1:m+1,1:m) = zeros(m+1,m);
cs(1:m) = zeros(m,1);
sn(1:m) = zeros(m,1);
e1    = zeros(n,1);
e1(1) = 1.0;


for iter = 1:max_it,                                    % begin iteration
     r = M \ ( b-A*x );   mvp_count = mvp_count+1;
     V(:,1) = r / norm( r );
     s = norm( r )*e1;
     for i = 1:m,                                       % construct orthonormal
         w = M \ (A*V(:,i));   mvp_count = mvp_count+1; % basis using Gram-Schmidt
         for k = 1:i,
             H(k,i)= w'*V(:,k);
             w = w - H(k,i)*V(:,k);
         end
         H(i+1,i) = norm( w );
         V(:,i+1) = w / H(i+1,i);
         for k = 1:i-1,                              % apply Givens rotation
             temp     =  cs(k)*H(k,i) + sn(k)*H(k+1,i);
             H(k+1,i) = -sn(k)*H(k,i) + cs(k)*H(k+1,i);
             H(k,i)   = temp;
         end
         [cs(i),sn(i)] = rotmat( H(i,i), H(i+1,i) ); % form i-th rotation matrix
         temp   = cs(i)*s(i);                        % approximate residual norm
         s(i+1) = -sn(i)*s(i);
         s(i)   = temp;
         H(i,i) = cs(i)*H(i,i) + sn(i)*H(i+1,i);
         H(i+1,i) = 0.0;
         error  = abs(s(i+1)) / bnrm2;
         if ( error <= tol ),                        % update approximation
            y = H(1:i,1:i) \ s(1:i);                 % and exit
            x = x + V(:,1:i)*y;
            break;
         end
     end

     if ( error <= tol ), break, end
     y = H(1:m,1:m) \ s(1:m);
     x = x + V(:,1:m)*y;                            % update approximation
     r = M \ ( b-A*x ); mvp_count = mvp_count+1;    % compute residual
     resvec(iter) = r;                              % store the residual AT the iter's iteration
     mvpvec(iter) = mvp_count;                      % store ACTUALLY performed MVP count BY iter's iteration
     s(i+1) = norm(r);
     error = s(i+1) / bnrm2;                        % check convergence
     if ( error <= tol ), break, end;
end

if ( error > tol ) flag = 1; end;                 % converged
% END of gmres.m