eig_lanczos.m

用lanczos求特征值的matlab程序

function [miu,x,qx,tol,k]=eig_lanczos(A)
k=0;
[w,e]=size(A);
r0=ones(w,1);
r0r0=norm(r0);
q1=r0/r0r0;
tol1(1)=1;
q0=0;
beta=1;
l=5;
y=0;
while tol1(k+1)>1e-4
    clear y;
    aq=A*q1;
    r=aq-beta*q0;
    alpha=q1'*r;
    r=r-alpha*q1;
    k=k+1;
    if k == 1
        T = [alpha];
        Q = [q1]; 
    elseif k == 2
        T = [T,beta;beta,alpha];
        Q = [Q,q1];
    else
        T1 = [zeros(k-2,1);beta];
        T = [T,T1;T1',alpha];
        Q = [Q,q1];
     end;
     tol1(k+1)=tol1(k);
        %if r~=0;
        beta=norm(r);
        q0=q1;
        q1=r/beta;
    %end
  if k>l
     h=1;
     T2=T(1:(k-1),:);
     T3=T2(:,1:(k-1));
     [v,d]=eig(T3);
     [m,n]=size(d);
     d=diag(d);
     for j=1:l
         for r=1:m
             if j==1
                d1=d;
             end
             a=max(d1);
             b=min(d1);
             if d(r)==a;
                y(:,h)=v(:,r);
                miu(k-l,h)=d(r);
                h=h+1;
                d1(r)=b;
             end
         end
     end
     Q1=Q(:,1:(k-1));
     y=y(:,1:l);
     miu=miu(:,1:l);
     for s=1:l
         x(:,s)=Q1*y(:,s);
     end
     beta=T(k,k-1);
     for t=1:l
         r(:,t)=beta*abs(y(m,t));
     end
     tol1(k+1)=min(abs(r));
     qt=Q1'*Q1-diag([ones(1,k-1)]);
     qx(k-l,1)=max(max(abs((qt))));
 end
end
k=k-l;
[g,e]=size(tol1);
tol=tol1(l+2:e);