function [A,flops]=QRitrid(A) % assume A is tridiagonal and symmetric n=size(A,1); i=n; % i=eigenvalue being computed normA=max(abs(diag(A)+[0; diag(A,1)]+[diag(A,1); 0])); % this is the infinity norm of A % takes 2n flops to compute flops=2*n; while i>1 sigma=A(i,i); % this is our shift B=A; [Q,R]=qr(B(1:i,1:i)-sigma*eye(i)); B(1:i,1:i)=R*Q+sigma*eye(i); % subtracting the shift from the diagonal of A (A=A-sigma*I) for j=1:i A(j,j)=A(j,j)-sigma; end for j=1:i-1 % figure out our Givens rotation x=A(j,j-(j>1)); y=A(j+1,j-(j>1)); c=x/sqrt(x^2+y^2); s=y/sqrt(x^2+y^2); % applying the Givens rotation on the left % on columns j-1 through j+2 st=max(1,j-1); en=min(j+2,i); A([j j+1],st:en)=[c s; -s c]*A([j j+1],st:en); % we multiply a 2x2 Givens rotation by a 2x4 matrix % we need to compute 8 elements at 3 flops per element % completing the similarity transformation on the right A(st:en,[j j+1])=A(st:en,[j j+1])*[c -s; s c]; end % adding back sigma * I for j=1:i A(j,j)=A(j,j)+sigma; end if abs(A(i,i-1))<10^(-16)*normA % eigenvalue #i is revealed so we deflate it i=i-1; disp('Moving on to eigenvalue number '); disp(i); end end z=diag(A);