% function a=SREigenvalues(B);
%
% Computes the eigenvalues of a sign regular matrix A
% A=CJ is assumed to be the product of a totally nonnegative matrix C
% and the reverse identity J(i,j)=delta_{i,n-j+1} 
% The input is the bidiagonal decomposition of C stored in B
%
% Examples:
% 1. SREigenvalues(ones(5)) computes the eigenvalues of the 
%    column-reversed 5-by-5 Pascal matrix
%
% 2. SREigenvalues(TNVandBD([0.1 0.2 0.4])) computes the eigenvalues 
%    of the row reversed Vandermonde matrix with nodes 0.1, 0.2, 0.4,
%    or (which is the same) the eigenvalues of the Vandermonde matrix
%    with nodes 0.4, 0.2, 0.1
%
% Reference: Plamen Koev and Froilan Dopico,
% Accurate Eigenvalues of Certain Sign Regular Matrices, Linear Algebra
% Appl.
%
% Copyright (c) 2006 Plamen Koev. See COPYRIGHT.TXT for more details.
% Written May 2006

function a=SREigenvalues(B);

n=size(B,1);

% reducing to upper triangular form
for i=1:n-1
  for j=n:-1:i+1
    x=B(j,i);
    B(j,i)=0;
    B=TNAddToNext(B',x,n-j+2)';
  end
end

% B is reduced to upper triangular form

% killing all but the diagonal and the superdiagonal

for i=n:-1:3
    for j=1:i-2
        % extracting B(j,i) all the way of the left
        x=B(j,i); 
        B(j,i)=0;
        for k=i-1:-1:j+1
            z=B(j,k)/(B(j+1,k+1)+x);
            B(j,k)=B(j+1,k+1)*z;
            B(j+1,k+1)=B(j+1,k+1)+x;
            x=z*x;
        end
        x=x*B(j,j)/B(j+1,j+1);
        B=TNAddToPrevious(B,x,1,n-j+1);
        
        % killing the below-diagonal buldge now
        x=B(n-j+1,n-j);
        B(n-j+1,n-j)=0;
        B=TNAddToNext(B',x,j+1)';
    end
end

% expanding the bidiaognal decomposition of a bidiagonal into an actual
% bidiagonal

for i=1:n-1
    B(i,i+1)=B(i,i+1)*B(i,i);
end

% symmetrizing the anti-bidiagonal matrix now

for i=1:fix(n/2)
    z=sqrt(B(i,i)*B(n-i+1,n-i+1));
    B(i,i)=z;
    B(n-i+1,n-i+1)=z;
    z=sqrt(B(i,i+1)*B(n-i,n-i+1));
    B(i,i+1)=z;
    B(n-i,n-i+1)=z;    
end

% the absolute values of the eigenvalues are the singular values of the
% column-reversed anti-bidiagonal

a=mexdlasq1(diag(B),diag(B,1));


% the signs of the eigenvalues are known from theory:

for i=2:2:n
    a(i)=-a(i);
end