How to reproduce the numerical examples in Prof. Saad's Book about Krylov subspace methods?

Question

After reading Prof. Saad' Book, "Iterative methods for Sparse Linear Systems, 2nd version", I want to do the numerical examples about the Krylov subspace methods not only to reproduce the results in his book but also to guarantee me really understand the idea of Krylov subspace methods, but I have some trouble.

First, in Page 98, section 3.7, he uses 5 matrices to test all his methods, the first 3 matrices are named F2DA,F2DB,F3D are from a software "SPARSKIT" using Fortran and under UNIX system, which is impossible for me to generate the 3 matrices, (I just can write matlab code under windows). so, I want to ask that whether anyone has done the numerical examples with these 3 matrices? Where can I download these 3 matrices? Or how to generate these 3 matrices?

Second, the last 2 matrices are from matrix market, named ORS and FID . So, I can only get these 2 matrices and use these 2 matrices to do numerical examples. (I feel that if an author uses the matrix from matrix market, it is so convient for us to get easily. ;). So I use the matrix ORS to do the first 2 experiments about full orthogonalization method (FOM) and GMRES on Page 168 and 180 respectively as follows:

Page 168 as follows: Example 6.1. Table 6.1 shows the results of applying FOM(10) with no preconditioning to three of the test problems described in Section 3.7. The column labeled Iters shows the total actual number of matrix-vector multiplications (matvecs) required to converge. The stopping criterion used is that the 2-norm of the residual be reduced by a factor of 107 relative to the 2-norm of the initial residual. A maximum of 300 matvecs are allowed. Kflops is the total number of floating point operations performed, in thousands. Residual and Error represent the two-norm of the residual and error vectors, respectively. Note that the method did not succeed in solving the third problem.

Table 6.1 A test run of FOM with no preconditioning \begin{array}{ccccc} \hline Matrix&Iters&Kflops&Residual&Error\\ F2DA&109&4442&0.36E-03&0.67E-04\\ F3D&66&11664&0.87E-03&0.35E-03\\ ORS&300&13558&0.26E+00&0.71E-04\\ \hline \end{array} Page 180, Table 6.2 A test run of GMRES with no preconditioning

Example 6.2. Table 6.2 shows the results of applying the GMRES algorithm with no preconditioning to three of the test problems described in Section 3.7. See Example 6.1 for the meaning of the column headers in the table. In this test, the dimension of the Krylov subspace is m = 10. Observe that the problem ORS, which could not be solved by FOM(10), is now solved in 205 steps.

\begin{array}{c|c|c|c|c} \hline Matrix&Iters&Kflops&Residual&Error\\ F2DA&95&3841&0.32E-02&0.11E-03\\ F3D&67&11862&0.37E-03&0.28E-03\\ ORS&205&9221&0.33E+00&0.68E-04\\ \hline \end{array}

My question is that when I write the FOM(10) method in matlab 2018b, I really cannot obtain the same resluts (Example 6.1) with that of his book (and I am sure my code is correct). In addition, I use the matlab built-in function gmres.m to test the result (example 6.2). I also cannot obtain the same results with his book. As we know, in FOM and GMRES method, one iteration step will require 1 matrix-vector, so the number of matrix-vector multiplication can be obtained by the iteration step. I want to know that how to obtain the results in his book i.e., What should I do to obtain the same results with his book? or indeed, I may not get the same results with his book if using matlab? or Indeed, I do not have to obtain the same results as long as my matlab code can converge? Thanks very much. below is my FOM(restart) matlab code.

% 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(restart) function is

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

Below is my experiments results:

************  FOM(restart) method :  ****************
Iters =
   300
Residual =
   6.9544e+03
Error =
  187.4011

************ gmres(restart)method:   ****************
Iters =
   300
Residual =
  216.7682
Error =
   24.0811