arnoldim_mine.m

来自「计算大型矩阵重启动的arnoldi的matlab 程序。」· M 代码 · 共 61 行

function [x,k,toltol] = arnoldim_mine(A,b)
[N,M] = size(A);
u =1e-8
m=20
aa = norm(A,1);
bb = norm(b);

x0 = 0;
r0 = b;
for k =1:50
    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,1)];
        if j == 1
            H = [h(j,j)];
            V = v(:,j);
        else
            H = [[H;zeros(1,j-2),h(j,j-1)],h(1:j,j);];
            V = [V,v(:,j)];
        end;
        %if rcond(H) < 1e16
        yk = H\rr;
        %    tol(j) = h(j+1,j)*abs(yk(j))/bb;
        %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;
    end;
   % x =x0 + V*yk;
    %tol(k) = norm(b - A*x)/(aa*norm(x)+bb);
    
    if tol((k-1)*m+j)< u
        disp('收敛退出');
        break;
    end;
    r0 = b - A*x(:,j);
    x0 =x(:,j);
end;
toltol = tol((k-1)*m+j);
k = (k-1)*m+j;
x=x(:,j);
save tol;
save x;