gmresm_mine.m

来自「计算大型矩阵的一般的gmres算法matlab程序。」· M 代码 · 共 60 行

function [x,k,tol1tol1] = gmresm_mine(A,b)
[M,N] = size(A);
u =1e-8
m=50

aa = norm(A,1);
bb = norm(b);

r0 = b;
x0 = 0;
for k = 1:100
r0r0 = norm(r0);
v(:,1) = r0/r0r0;
H =zeros(0,0);
for j =1:m
    w = A*v(:,j);
    for i =1:j
        h(i,j) = v(:,i)'*w;
        w = w - h(i,j)*v(:,i);
    end;
    h(j+1,j) = norm(w);
    v(:,j+1) = w/h(j+1,j);
    rr = r0r0*[1;zeros(j,1)];
    if j == 1
        H = [h(j,j);h(j+1,j)];
        V = v(:,j);
    else
        H = [H,h(1:j,j);zeros(1,j-1),h(j+1,j)];
        V = [V,v(:,j)];
    end;
    
    [Q1,R] = qr(H);
    gk = Q1'*rr;
    R = R(1:j,:);
    gk = gk(1:j,:);
    yk = R\gk;
    %if mod(j,10)==0
    %x(:,j) = x0 + V*yk;
    %tol((k-1)*m+j) = norm(b - A*x(:,j))/(aa*norm(x(:,j))+bb);
    %if tol((k-1)*m+j)< u
        %break;
        %end;
    %end;
end;
x = x0 + V*yk;
tol(k) = norm(b - A*x)/(aa*norm(x)+bb);
if tol(k)< u
        disp('收敛退出');
        break;
end;
r0 = b - A*x;
x0 =x;
end;
tol1tol1 = tol(k);
save tol;