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
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. $\endgroup$ – Lutz Lehmann Dec 14 '19 at 12:53