%Procedure RSToepeig.m
%Copyright 1997 William F. Trench
%
%
%As defined by A. L. Andrew, a vector $x=[x_{0} x_{1} \cdots x_{n-1}]$
%is said to be symmetric (skewsymmetric) if $x_{r}=x_{n-r}$
%($x_{r}=-x_{n-r}$), $0 \le r \le n-1$. A $n=2m$ a real symmetric
%Toeplitz matrix $T=[t_{|r-s|}]_{r,s=0}^{n-1}$ has $m$ linearly
%independent symmetric eigenvectors and $m$ linearly independent
%skewsymmetric eigenvectors. A $n=2m+1$ then $T$ has $m+1$ linearly
%independent symmetric eigenvectors and $m$ linearly independent
%skewsymmetric eigenvectors. As defined by Trench, an eigenvalue is
%"even" (parity $0$) if it is associated with a symmetric eigenvector
%or odd (parity $1$) if it is associated with a skew symmetric
%eigenvector. By a theorem of Genin, if $\lambda$ is an eigenvalue of
%multiplicity $2r$ then the $\lambda$-eigenspace of $T$ has $r$
%linearly independent symmetric eigenvectors and $r$ linearly
%independent skew symmetric eigenvectors, and if $\lambda$ has
%multiplicity $2r+1$ then the $\lambda$-eigenspace of $T$ has a basis
%consisting of either $r+1$ symmetric and $r$ skew symmetric
%eigenvectors, or of $r$ symmetric and $r+1$ skewsymmetric
%eigenvectors. Thus, a repeated eigenvalue is necessarily both even
%and odd.
%
%This procedure computes the eigenvalues of a real symmetric Toeplitz
%matrix with first column vector $tvec$. The even eigenvalues are the
%elements of the vector $Mu$ and the associated normalized eigenvectors
%are the columns of $P$. The odd eigenvalues are the elements of the
%vector $Nu$ and the associated normalized eigenvectors are the columns
%of $Q$. In addition, $L = [Mu' Nu']'$ and $Omega = [P Q] are
%available. The eigenvalues are presented in nondecreasing order with
%with their parities (parity = 0 for even, 1 for odd).
%
%To execute the procedure, first enter $tvec$ by declaring
% "tvec= (column vector) Then enter "RSToepeig" (not case-sensitive).
%
%To compute the Jacobian of the eigenvalues with respect to
%$t_0,...t_(n-1)$, define "compjac=[]" prior to starting.

 [n,x] = size(tvec);
 r = ceil(n/2);
 s = floor(n/2);
 A = zeros(r,r);
 B = zeros(s,s);
 % If the order of $T$ is even go this way.
 if r == s
     for i = 1:r
           for j = i:r
                A(i,j) = tvec(j-i+1);
                B(i,j) = tvec(i+j);
           end
           for j = 1:i-1
                A(i,j) = A(j,i);
                B(i,j) = B(j,i);
           end
     end

     [P0,E] = eig(A+B);       %even eigenvalues
     P = ((2)^(-1/2))*[(flipud(P0))',P0']';  %corresponding eigenvectors
     [Q0,O] = eig(A-B);       %odd eigenvectors
     Q = (2^(-1/2))*[-(flipud(Q0))',Q0']';   %corresponding eigenvectors
     evenoddeigs =  [sort(diag(E))' sort(diag(O))']';



 else
 % If the order of $T$ is odd go this way.

     A = zeros(r,r);
     A(1,1) = tvec(1);
     for j = 2:r
           A(1,j) = 2^(1/2)*tvec(j);
     end
     for i = 2:r
           for j = i:r
                A(i,j) = tvec(j-i+1)+tvec(j+i-1);
           end
           for j = 1:i-1
                A(i,j) = A(j,i);
           end
     end


     [P0,E] = eig(A);  %even eigenvalues
     P1=P0;
     P1([1],:)=[];
     P = (2^(-1/2))*[(flipud(P1))',2^(1/2)*P0([1],:)',P1']';
                    %corresponding eigenvectors

     for i = 1:s
           for j = i:s
                B(i,j) = tvec(j-i+1)-tvec(j+i+1);
           end
           for j = 1:i-1
                B(i,j) = B(j,i);
           end
     end


     [Q0,O] = eig(B);    %odd eigenvalues
     Q = (2^(-1/2))*[-(flipud(Q0))',zeros(s,1),Q0']';
                    %corresponding eigenvectors
 end

 [Mu, order] = sort(diag(E));
 P=P(:,order);
 [Nu, order] = sort(diag(O));
 Q=Q(:,order);

  eigenvalues = [Mu' Nu']';
  L = [Mu' Nu']';

  Omega = [P Q];

  parity= [zeros(1,ceil(n/2))  ones(1,floor(n/2))]';
  eigenvalues1 = eigenvalues;
  for i = 1:n-1
       for j = i+1:n
          if eigenvalues1(j) < eigenvalues1(i)
             tempeig=eigenvalues1(i);
             temppar=parity(i);
             eigenvalues1(i)=eigenvalues1(j);
             parity(i)=parity(j);
             eigenvalues1(j)=tempeig;
             parity(j)=temppar;
          end
       end
  end

  eigenvalues = sort(eigenvalues);

  if exist('noprint') == 0
  eigenvalues = sort(eigenvalues);
  fprintf('\n                eigenvalue      parity\n\n')
  for i = 1:n fprintf('          %+15.14e   %g\n',eigenvalues(i), parity(i)),end
  fprintf('\n')
  end

%Compute the jacobian of the eigenvalues w.r.t. $t_0,...,t_(n-1)$


 if  exist('compjac') == 1
     C = zeros(n,n);
     for i=1:n   C(i,1)=1; C(i,n) =2*Omega(1,i)*Omega(n,i); end
     for j=2:n-1
       mm = ceil((n-j)/2);
       mmm = floor((n-j)/2);
          if mm == mmm
                  for i=1:n
                      C(i,j) = 4*Omega(j:j+mm-1,[i])'*Omega(1:mm,[i])+2*Omega(j+mm,i)*Omega(mm+1,i);
                  end
         else
                  for i = 1:n
                      C(i,j) = 4*Omega(j:j+mm-1,[i])'*Omega(1:mm,[i]);
                 end
        end
     end
     fprintf(' Condition number of Jacobian matrix =  %3.2e  \n', cond(C))
   end
