在阅读了 Saad 教授的书“稀疏线性系统的迭代方法，第 2 版”之后，我想做关于 Krylov 子空间方法的数值示例，不仅可以重现他书中的结果，而且可以保证我真正理解这个想法克雷洛夫子空间方法，但我有一些麻烦。
首先，在第 98 页，第 3.7 节中，他使用 5 个矩阵来测试他的所有方法，前 3 个矩阵名为F2DA，F2DB，F3D来自使用 Fortran 和 UNIX 系统的软件“ SPARSKIT ”，这对我来说是不可能的生成3个矩阵，（我只能在windows下编写matlab代码）。所以，我想问一下是否有人用这三个矩阵做过数值例子？我在哪里可以下载这 3 个矩阵？或者如何生成这3个矩阵？
其次，最后两个矩阵来自矩阵市场，分别命名为ORS和FID。所以，我只能得到这2个矩阵，并用这2个矩阵来做数值例子。（我觉得如果作者使用矩阵市场的矩阵，我们很容易得到。;）。所以我使用矩阵ORS分别在第 168 页和第 180 页上进行了关于全正交化方法 (FOM) 和 GMRES的前 2 个实验，如下所示：
第168页如下：示例6.1。表 6.1 显示了将 FOM(10) 应用到第 3.7 节中描述的三个测试问题的结果。标记为 Iters 的列显示了收敛所需的矩阵向量乘法 (matvecs) 的实际总数。使用的停止标准是残差的 2 范数相对于初始残差的 2 范数减少了 107 倍。最多允许 300 个 matvecs。Kflops 是执行的浮点运算的总数，以千计。残差和误差分别表示残差和误差向量的二范数。请注意，该方法未能成功解决第三个问题。
表 6.1 无预处理的 FOM 测试运行
第 180 页，表 6.2 没有预处理的 GMRES 测试运行
例 6.2。表 6.2 显示了对 3.7 节中描述的三个测试问题应用不带预处理的 GMRES 算法的结果。有关表中列标题的含义，请参见示例 6.1。在这个测试中，Krylov 子空间的维数是 m = 10。观察到，FOM(10) 无法解决的 ORS 问题现在分 205 步解决。
我的问题是，当我在 matlab 2018b 中编写 FOM(10) 方法时，我真的无法获得与他的书相同的结果（示例 6.1）（我确信我的代码是正确的）。另外，我使用了matlab内置函数gmres.m测试结果（示例 6.2）。我也无法用他的书获得相同的结果。众所周知，在 FOM 和 GMRES 方法中，一个迭代步需要 1 个矩阵向量，因此矩阵向量乘法的次数可以通过迭代步得到。我想知道如何获得他书中的结果，即，我应该怎么做才能获得与他的书相同的结果？或者实际上，如果使用 matlab，我可能不会得到与他的书相同的结果？或者确实，只要我的matlab代码可以收敛，我就不必获得相同的结果？非常感谢。下面是我的 FOM（重启）matlab 代码。
% main.m
clc;clear;close all;
% load the sparse matrix 'ORS' from matrix market and convert it to matlab format
load orsirr_1.mtx
i = orsirr_1(2:end,1);
j = orsirr_1(2:end,2);
elem= orsirr_1(2:end,3);
A = sparse(i,j,elem);
% initizlization
n = length(A);
x_star = ones(n,1);% set the exact solution is [1,...,1]'
b = A*x_star;% generate the right hand side of Ax=b
tol = 1e-7;
maxit = 30;% maximum the outer iteration
restart = 10;% restarted steps
x0 = zeros(n,1);% initial guess vector
%% FOM(restart)
fprintf('\n************ FOM(restart) method : ****************\n')
[x,flag,relres,iter,resvec] = myFOM_restart(A,b,tol,maxit,x0,restart);
Iters=iter, Residual =resvec(end), Error = norm(x_star-x)
%% gmres(restart)
fprintf('\n************ gmres(restart)method: ****************\n')
[x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,[],[],x0);
Iters=(iter(1)-1)*restart+iter(2);
Iters, Residual =resvec(end), Error = norm(x_star-x)
FOM（重启）功能是
function [x,flag,relres,iter,resvec] = myFOM_restart(A,b,tol,maxit,x0,restart)
% Full Orthogonalization Method (FOM(m)) using modified Gram-Schmidt with
% restarted solving Ax=b
% input
% A real nonsingular matrix in n X n
% b right hand side
% tol tolerance of relative residual norm: ||r_m||/||r_0||<tol
% maxit outer maxmimu iterations
% restart restarted steps
% x0 initial guess vector usually taken zero
% output
% x approximate solution
% flag convegence :0
% unconvergence 1
% relres relative residual norm ||r_m||/||r_0||
% iter iteration steps
% resvec all the residual vector norm [||r0||,...||rm||]
% written on 2019 12.14
% reference P167 Algorithm 6.4.1 from
% Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia,
% PA, second edition, 2003.
%% initialize space to guarantee efficience
m = restart;
n = length(b);
H = zeros(m+1,m);
resvec = zeros(m+1,1);
V = zeros(n,m+1);
r0 = b-A*x0;
beta = norm(r0);
V(:,1) = r0/beta;
resvec(1) = beta;
flag = 1;
%% begin to iteration
for k=1:maxit% loop of maximum outer iteration steps
% modified Gram-Schmidt orthogonalization
for j = 1:m
w = A*V(:,j);
for i=1:j
H(i,j) = w'*V(:,i);
w = w-H(i,j)*V(:,i);
end
H(j+1,j) = norm(w);
if H(j+1,j) ==0
flag = 0;
fprintf('lucky breakdown!!!!!!!\n')
break;
end
V(:,j+1) = w/H(j+1,j);
% compute the new residual vector
e1 = zeros(j,1);e1(1)=1;
ej = zeros(j,1);ej(j)=1;
y = H(1:j,1:j)\(beta*e1);% ym
resvec(j+1,1) = norm(-H(j+1,j)*ej'*y*V(:,j+1));% norm of residual r_j
% check convergence
relres = resvec(j+1,1)/norm(r0);
if relres < tol
flag = 0;
break;%
end
end%end of restart
if flag==0
break;
end
% update new approximate solution
x = x0+V(:,1:j)*y;
x0 = x;
end % end of outer maximum iteration
iter = j+(k-1)*restart;% total iteration steps
if flag==0
x = x0+V(:,1:j)*y;
fprintf('** FOM(m) converges in %d steps to relative residual %5.2e*****\n',iter,relres);
else
% fprintf('** FOM(m) unconverges in %d steps with relative residual %5.2e *\n',iter,relres);
end
resvec(iter+2:end) = [];
end
下面是我的实验结果：
************ FOM(restart) method : ****************
Iters =
300
Residual =
6.9544e+03
Error =
187.4011
************ gmres(restart)method: ****************
Iters =
300
Residual =
216.7682
Error =
24.0811