function [a,b] = qrStep(a,b)
% Computes the orthogonal similarity tranformation Q^TAQ for a symmetric
% tridiagonal matrix A, using rotation matrices. 
% Here Q is the orthogonal matrix from the QR-factorization of A.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
% a: The main diagonal of A. Vector of length n.
% b: The sub/superdiagonal of A. Vector of length n-1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT
% a: The main diagonal of Q^TAQ. Vector of length n.
% b: The sub/superdiagonal of Q^TAQ. Vector of length n-1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Check lengths
n = length(a);
if length(b) ~= n-1
    error('The sub/superdiagonal must have one less element than the diagonal')
end


% Allocate space for the sine and cosine of the rotation angle
c = zeros(n-1,1);
s = zeros(n-1,1);

% Initialize a temporary variable 
t = b(1);

% Compute (part of) the QR factorization
for j = 1:n-1
    
    % Compute the required roation angle to make element (j+1,j) of A zero.
    r = sqrt(a(j)^2+t^2);
    c(j) = a(j)/r;
    s(j) = t/r;
    
    % Multiply implicitly by the corresponding rotation matrix.
    a(j) = r;
    t = b(j);
    b(j) = t*c(j)+a(j+1)*s(j);
    a(j+1) = -t*s(j)+a(j+1)*c(j);
    if j < n-1
        t = b(j+1);
        b(j+1) = t*c(j);
    end
    
end


% Compute the non-zero elements of RQ = Q^TAQ
for j=1:n-1
    
    % Multiply with the transpose of the j-th rotation matrix computed
    % earlier
    a(j) = a(j)*c(j)+b(j)*s(j);
    b(j) = a(j+1)*s(j);
    a(j+1) = a(j+1)*c(j);
end

end