chenmh_arnoldi.m

来自「Frequently used algorithms for numerical」· M 代码 · 共 33 行

function [x,k,err] = Chenmh_Arnoldi(A,b,x0,tol,mode);
m = size(A,1);
r0 = b-A*x0; nr0 = (r0'*r0)^0.5;
err = [];
if (nr0==0)
    x = x0;
else
    V = r0/nr0; H = []; k = 0; istop = 0;
    while ( (k<m) & (istop==0) )
        [V,H,k,istop] = Chenmh_ArnoldiPro(A,m,k,V,H);
        e1 = eye(k); e1 = e1(:,1);
        if (mode==0) % Arnoldi
            y = H(1:end-1,:)\e1*nr0;
            err(k) = H(end,end)*abs(y(end));
        else         % GMRES   
            [Q,R] = qr(H);                                
            g = Q'*nr0*[e1;0];
            y = R(1:end-1,:)\g(1:end-1);
            err(k) = abs(g(end));
        end
        if (err(k)<tol)
            istop = 1;
        end
    end
    x = x0+V(:,1:end-1)*y;
end