GMRES:
function [ x,j,res,resvec,time] = mygmres( A, b, x0, tol, max_it )
%% GMRES METHOD
% j : the number of iterations;
% res : the residual value at termination;
% x : solution
% time : CPU time;
if nargin < 5
max_it=1000;
end
if nargin < 4
tol = 1.e-10;
end
if nargin < 3
x0 = zeros(length(b),1);
end
tic;
r0 = A*x0 - b;
beta = norm(r0);
V(:,1) = r0/beta;
resvec = beta;
Q = cast(1,'like',b);
j = 0;
while ((abs(resvec(j + 1)) > tol) && (j < max_it))
% Arnoldi %
j = j + 1;
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,2);
V(:,(j + 1)) = w/H((j + 1),j);
% Construct R and Givens rotation %
H(1:j,j) = Q*H(1:j,j);
rho = H(j,j);
H(j,j) = sqrt(rho^2 + H((j +1),j)^2);
c = rho/H(j,j);
s = H((j + 1),j)/H(j,j);
H((j + 1),j) = 0;
% Apply Givens rotation to Q %
Q((j + 1),:) = -s*Q(j,:);
Q(j,:) = c*Q(j,:);
Q((j + 1),(j + 1)) = c;
Q(j,(j + 1)) = s;
% Apply Givens rotation to resvec%
resvec(j + 1,1) = -s*resvec(j,1);
resvec(j,1) = c*resvec(j,1);
end
y = H((1:j),:)\resvec(1:j);
x = x0 - V(:,(1:j))*y;
res = abs(resvec(j + 1));
resvec = abs(resvec);
time=toc;
end
A=[ 20, -10., 0., 0., 0., 0., 0., 0., 0.;
-10., 20., -10., 0., 0., 0., 0., 0., 0.;
0., -10., 20., -10., 0., 0., 0., 0., 0.;
0., 0., -10., 20., -10., 0., 0., 0., 0.;
0. , 0., 0., -10., 20., -10., 0., 0., 0.;
0., 0., 0., 0., -10., 20., -10., 0., 0.;
0., 0., 0., 0., 0., -10., 20., -10., 0.;
0., 0., 0., 0., 0., 0., -10., 20, -10.;
0., 0., 0., 0., 0., 0., 0., -10., 20.]
b=[0.11407755, 0.048673, 0.03304509, 0.02532112, 0.02064553, 0.01749003, 0.01520812, 0.01347669, 0.0121155 ]'
>> mygmres( A, b, x0, 2*10^(-8), 9 )
ans =
0.0206
0.0297
0.0340
0.0350
0.0335
0.0299
0.0245
0.0177
0.0094
>> b'/A
ans =
0.0206 0.0297 0.0340 0.0350 0.0335 0.0299 0.0245 0.0177 0.0094
预优GMRES:
function [ x,j,res,resvec,time] = mygmres_ja( A, b, x0, tol, max_it )
%% GMRES METHOD (Pre-processing M = Jacobi)
% j : the number of iterations;
% res : the residual value at termination;
% x : solution
% time : CPU time;
if nargin < 5
max_it=1000;
end
if nargin < 4
tol = 1.e-10;
end
if nargin < 3
x0 = zeros(length(b),1);
end
tic;
M = diag(diag(A)); % Jacobi pre-processing matrix;
r0 =M\(A*x0 - b);
beta = norm(r0);
V(:,1) = r0/beta;
resvec = beta;
Q = cast(1,'like',b);
j = 0;
while ((abs(resvec(j + 1)) > tol) && (j < max_it))
% Arnoldi %
j = j + 1;
w = M\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,2);
V(:,(j + 1)) = w/H((j + 1),j);
% Construct R and Givens rotation %
H(1:j,j) = Q*H(1:j,j);
rho = H(j,j);
H(j,j) = sqrt(rho^2 + H((j +1),j)^2);
c = rho/H(j,j);
s = H((j + 1),j)/H(j,j);
H((j + 1),j) = 0;
% Apply Givens rotation to Q %
Q((j + 1),:) = -s*Q(j,:);
Q(j,:) = c*Q(j,:);
Q((j + 1),(j + 1)) = c;
Q(j,(j + 1)) = s;
% Apply Givens rotation to resvec%
resvec(j + 1,1) = -s*resvec(j,1);
resvec(j,1) = c*resvec(j,1);
end
y = H((1:j),:)\resvec(1:j);
x = x0 - V(:,(1:j))*y;
res = abs(resvec(j + 1));
resvec = abs(resvec);
time=toc;
end
见《https://blog.csdn.net/waitingwinter/article/details/103863358》