% 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