从技术上讲， Netlib提供的模板 Matlab 代码已经计算了每次迭代获得的残差。它只是没有以这是子例程的输出的方式记录。
每次迭代后都会计算iter残差r，因此只需记录一下即可。在残差向量中计算后我会立即执行此操作resvec(iter)。
现在，提供的 GMRES 代码不是最优的（例如，残差r是在循环之前和之后计算的 -i变量：正交基的构造），所以我不打算优化该行为，而只是计算矩阵向量产品 (MVP) 每次发生在mvp_count. 然后，在每次迭代结束时，我会将实际执行的 MVP 的数量记录在mvpvec(iter).
两者resvec和mvpvec都添加到gmres函数的输出参数中，该参数应提供足够的信息来进行绘图。
笔记：
mvpvec(iter)计算的特定实现在第 th 次迭代结束时执行的 MVP 。itergmres- 没有通过 RHS 的 2 范数进行归一化
b（以避免残差、相对误差等之间的混淆）。应该在函数内部或外部轻松更改。
- 只有具有所谓大原始矩阵的 MVPA正在被计算
- 不记录初始残差和相应的 MVP（在迭代过程开始之前）。
这是稍微修改的Matlab代码。
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