lanczos_mine.m

对称矩阵的lanczos算法matlab程序

function [x,j,toltol] = lanczos_mine(A,b)
[N,M] = size(A);
aa = norm(A,1);
x = zeros(N,1);
u =1e-18

%x = [1:1000]
%xx = norm(x)
bb = norm(b);
%aa = norm(A,1)
%axb = aa*xx+bb


%pause;
r0 = b;
r0r0 = norm(r0);

qj_1 = 0;
qj = r0/r0r0;
hj_1 = 1;

for j =1:250
    rj =A*qj-hj_1*qj_1;
    aj = qj'*rj;
    rj = rj - aj*qj;
    rr = r0r0*[1;zeros(j-1,1)];
    if j == 1
        T = [aj];
        Q = [qj]; 
    elseif j == 2
        T = [T,hj_1;hj_1,aj];
        Q = [Q,qj];
    else
        T1 = [zeros(j-2,1);hj_1];
        T = [T,T1;T1',aj];
        Q = [Q,qj];
    end;
    %if norm(rj)/bb > u
        hj = norm(rj);
        qj1 = rj/hj;
        if rcond(T) < 1e16
            yk = T\rr;
            x(:,j) =  Q*yk;
            tol(j) = norm(b - A*x(:,j))/(aa*norm(x(:,j))+bb);
            %tol(j) = hj*abs(yk(j))/bb;
            if tol(j)< u
                disp('收敛退出');
                break;
            end;
        end;
        %else
        %yk = T\rr;
        %disp('过程中断退出');
        %break;
        %end;
    qj_1 = qj;
    qj = qj1;
    hj_1 = hj;
end;
%i =[1:20]
%subplot(2,1,1)
%plot(i,i)
%subplot(2,1,2)
%xx = x(:,j);
toltol = tol(j);
save tol;
%save x;
%semilogy(i,tol(i))
x =  Q*yk;
T(200:240,200:240)