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

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

• What is the specific problem that you think requires UNIX? Nothing on the linked web page indicates that. The gzipped tape archive is not dependent on the operating system. The availability of FORTRAN compilers is also not restricted to specific OS. – Lutz Lehmann Dec 14 '19 at 10:53
• Thanks for your reply, I do not know about Fortran, I thought the SPARSKIT software can only work in UNIX. Though I have downloaded the software, I still do not know how to operate the Fortran software. Now, I just use the last 2 matrices to do numerical examples. Do you know how to reproduce the results with the author's results for ORS matric? Thanks. – sunshine Dec 14 '19 at 12:17
• I partially retract my comment, while there are Fortran compilers for MS Visual Studio etc., they are probably too expensive for just this purpose. To get the gnu toolchain with make etc. to compile this code, you will need a gnu environment, like cygwin or MinGW provide (no Linux or Unix kernel), or, as I said, some Linux virtual machine that can be minimal (gui-less or with some minimal windowing system like a tabbing ui or xfce) and thus quite small. // lfortran.org might become interesting in some distant future, "pre-alpha" hints at unexpected surprises. – Lutz Lehmann Dec 14 '19 at 12:53
• Thanks for your detailed explanation, I will try later when I complete the results with the last 2 matrices. Because My problem is now, I have got the ORS matrix and implement the example for FOM and gmres method. But I cannot get the results with the author's book. Because I think when I can reprodue the results of ORS matrix, it means that my matlab code is correct and it will work also for other matrics. So, I am firstly to debug my code for ORS matrix but I cannot find what is wrong with my code. Thanks very much. – sunshine Dec 14 '19 at 13:01