qr2shift.m

来自「various matlab programs to slove various」· M 代码 · 共 26 行

function [T,iter]=qr2shift(A,toll,itmax)
n=max(size(A)); iter=0; [T,Q]=houshess(A);
for k=n:-1:2
  I=eye(k); 
  while abs(T(k,k-1)) > toll*(abs(T(k,k))+abs(T(k-1,k-1)))
     iter=iter+1;  if (iter > itmax), return, end
     mu=T(k,k); [Q,R,c,s]=qrgivens(T(1:k,1:k)-mu*I);
     T(1:k,1:k)=R*Q+mu*I; 
     if (k > 2), 
        Tdiag2=abs(T(k-1,k-1))+abs(T(k-2,k-2));
        if abs(T(k-1,k-2)) <= toll*Tdiag2;
           [lambda]=eig(T(k-1:k,k-1:k));
           [Q,R,c,s]=qrgivens(T(1:k,1:k)-lambda(1)*I);
           T(1:k,1:k)=R*Q+lambda(1)*I;
           [Q,R,c,s]=qrgivens(T(1:k,1:k)-lambda(2)*I);
           T(1:k,1:k)=R*Q+lambda(2)*I;
        end
     end
  end, T(k,k-1)=0;
end
I=eye(2);
while (abs(T(2,1)) > toll*(abs(T(2,2))+abs(T(1,1)))) & (iter <= itmax)
   iter=iter+1; mu=T(2,2);
   [Q,R,c,s]=qrgivens(T(1:2,1:2)-mu*I);  T(1:2,1:2)=R*Q+mu*I;
end
return