Trench Real Symmetric Inverse Eigenvalue Solver.

This package consists of four m-files, separated by borders
like this.

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

They should be separated into separate files titled
1. ties.m 2.
2. tiessub.m
3. newties.m
4. toepinit.m

I make no claims for programming elegance in these  files
(except for toepinit.m, which was created by Dirk P. Laurie
and provides an elegant and efficient way to compute starting
vectors for iteration performed by ties.m); however, it works!

Contact me (wtrench@trinity.edu) if you have questions
concerning  thispackage.



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx




%Program ties.m (Toeplitz Inverse Eigenvalue Solver).
%William F. Trench
%Mathematics Department
%Trinity University
%San Antonio, Texas 78212 USA
%wtrench@trinity.edu

%This program attempts to find a vector toepgen which defines a
%Toeplitz matrix T with a given  spectrum target0. In the  default
%mode, T has  the same parity sequence as T_0, the Toeplitz
%matrix generated by the starting generator.

%It calls tiessub.m as a subroutine.

%tvec0 = generator for the initial n x n Toeplitz matrix T_0 (input).
%        The default generator comes from f(\theta)=-\theta
%target0 = desired spectrum (input) (column vector).
%parity0 = parity sequence of the eigenvalues of T_0 (computed below).

%The output is appended to the file ties.out.

%REF :=  The paper  "Numerical solution  of the  inverse eigenvalue problem
%for real symmetric Toeplitz matrices."

%Before starting  the program  it is  necessary to  define target0,
%the target spectrum, and  tvec0, the  generator of  the initial
%Toeplitz matrix  T_0. If tvec0 is not specified, the default generator
%is as defined above. (See  REF.) The target spectrum is normalized before
%the computation starts. The  solution of the normalized problem  is "tvec",
%which generates  an RST with trace  zero and Frobenius norm one.
%The solution to the original (unnormalized) problem is "toepgen".

%Phase 1 of the  computation is the determination  of the spectrum of
%T_0, the RST generated by tvec0. If you wish to use the same starting
%generator  to
%solve the IEP for several different  target spectra, then enter
%"repeat =  []" after the first  computation. To end  the "repeat" mode
%enter "clear repeat"; however, this is done automatically if you
%change n (order of the matrix).

%In  the  default  (automatic)  mode  the  iteration starts with rho=autorho
%(defined below as 0) and  linquit=autolinquit (defined below as 4).  If you
%wish  to  be  prompted  for  these  inputs,  then enter "enter = []" before
%starting  the  computation.  To  resume  the  automatic  mode, enter "clear
%enter".

%To start the program enter "ties".

%The final results are printed to "ties.out" in the current directory.
%If  this file does not exist then it is created. If it exists the results
%are  appended to it. In the default mode, only  the iteration
%count is printed to the  file. To  get  "tvec"  (the  answer  to
%the  normalized  problem)  printed,   enter "fileprint=[]" before starting
%the computation.

%Note: "toepgen" (the answer you probably want) is not printed to ties.out;
%to record toepgen, you can use "diary on" and enter "toepgen."
%To have  the final  results printed  to the  screen, enter
%"screenpint =  []" before starting the program. However, there  is little
%point to this if  n>20, since the results will simply scroll off the screen.
%To  plot  the   trigonometric  polynomial
%                  t_0+2\sum_{m=1}^{n-1}t_m\cos m(\pi*x)
%automatically at  the end  of the  computation, enter  "plotpoly =  []"
%before starting. If you decide to plot this polynomial after completion of
%ties, then enter "plotonly = []" and restart ties.





 if exist('plotonly')==0
 n = length(target0);
 r = ceil(n/2);
 s = floor(n/2);

if exist('tvec0') == 0
   clear repeat
   tvec0 = [];
end
if length(tvec0) ~= n clear repeat; end



  if exist('endquit')==0 endquit = 10^(-14); end %used to set accuracy
         %requirement for final answer
  if exist('jump1')==0 jump1 = 1;  end  % or 10^(-jump1)*sigma if rho > 0
  if exist('badmax')==0 badmax = 1; end
  if exist('deltarho')==0 deltarho = .1; end %automatic increment in rho
           %after failure with a given rho

  %rho and linquit are defined automatically if the variable 'enter' is undefined

  if exist ('autorho')==0 autorho = 0; end  %initial value of rho in automatic mode
  rho =autorho;
  if exist('autolinquit')==0 autolinquit = 4; end %value of linquit in automatic mode

  if exist('enter') == 1
       rho = input(' rho (0 <= rho < 1) = ');
       if isempty(rho) rho = autorho; end
  else
       rho = autorho;
  end

  if exist('enter') == 1
     fprintf('\n  Set rho = 0 when ||spectrum - target|| < 10^(-I) with I < 13\n')
     linquit = input(' I = ');
     if isempty(linquit) linquit = autolinquit; end
  else
     linquit = autolinquit;
  end
  if exist('fileprint') == 1 fprintf('ties.out','  n = %g   rho = %3.2f  linquit = %f\n',n, rho,linquit),end

  % default tvec0

 if exist('repeat') == 0

  fprintf('\n  Computing eigenvalues and eigenvectors of toep(tvec0).\n')


  if exist('tvec0') == 0
      tvec0=zeros(n,1);
      for i = 1:n
      if ceil(i/2) == floor(i/2) tvec0(i) = 1/(i-1)^2;end
      end
  end
  if length(tvec0) ~= n
      tvec0=zeros(n,1);
      for i = 1:n
        if ceil(i/2) == floor(i/2) tvec0(i) = 1/(i-1)^2;end
      end
  end

 A = zeros(r,r);
 B = zeros(s,s);

 % If n is even go this way.

 if r == s
     for i = 1:r
           for j = i:r A(i,j) = tvec0(j-i+1); B(i,j) = tvec0(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 = [(flipud(P0))',P0']';  %corresponding eigenvectors
     [Q0,O] = eig(A-B);       %odd eigenvectors
     Q = [-(flipud(Q0))',Q0']';   %corresponding eigenvectors
 else

 % If n is odd go this way.

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

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

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

     [P0,O] = eig(B);  clear B  %odd eigenvalues
     Q = [-(flipud(P0))',zeros(s,1),P0']';
                    %corresponding eigenvectors
 end

 [Even, order] = sort(diag(E));
 P=P(:,order);
 [Odd, order] = sort(diag(O));
 Q=Q(:,order); clear order

 %This finishes the computation of the eigenvectors.

 spectrum0 = [Even',Odd']';
%column vector with even eigenvalues in nondecreasing order, followed by odd
%eigenvalues in nondecreasing order.

 E1=diag(E); O1=diag(O);
 clear Even Odd  P0 Q0  E O

%Determine the parity sequence  of T_1:

  parity0= [zeros(1,ceil(n/2))  ones(1,floor(n/2))]';
  spectrum1 = spectrum0;
  for i = 1:n-1
       for j = i+1:n
          if spectrum1(j) < spectrum1(i)
             tempeig=spectrum1(i);
             temppar=parity0(i);
             spectrum1(i)=spectrum1(j);
             parity0(i)=parity0(j);
             spectrum1(j)=tempeig;
             parity0(j)=temppar;
          end
       end
  end
  clear spectrum1 tempeig temppar

     R = [P, Q];  clear P Q
% matrix of eigenvectors of T, each with norm=\sqrt2
     C0 = zeros(n,n);
     for i=1:n   C0(i,1)=1; C0(i,n) = R(1,i)*R(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
                      C0(i,j) = 2*R(j:j+mm-1,[i])'*R(1:mm,[i])+R(j+mm,i)*R(mm+1,i);
                  end
         else
                  for i = 1:n
                      C0(i,j) = 2*R(j:j+mm-1,[i])'*R(1:mm,[i]);
                 end
        end
     end
 end % of if exist('repeat') == 0

  delta = 2500; %initialize to enable restart if terminal computation
      %does not converge after \rho is reduced to zero
  rhoswitch = 0; %initialize switch to save unnecessary computation of
          %sigma when reducing \rho to zero for terminal computation
  rholast = 0;
  badcount = 0;
  tvec = tvec0;
  parity = parity0;
  spectrum =  spectrum0; % input for first pass through tiessub.m

%Develop  target  by  normalizing  target0;  target  will have the following
%properties:

%(a) The sum of its elements = 0.
%(b) Its norm = 1.
%(c) It has the same parity sequence as the eigenvalues of T_1

  target  = target0 - (sum(target0)/n)*ones(n,1);
  target = target/norm(target);

 %order target so that the eigenvalues of T have the same parity
 %sequence as the eigenvalues of T_0
  temptarget = sort(target);
  Even = zeros(1,r);
  Odd = zeros(1,s);
  Evencount = 1;
  Oddcount = 1;
  for i = 1:n
      if parity0(i) == 0
          Even(Evencount) = temptarget(i);
          Evencount = Evencount + 1;
      else
          Odd(Oddcount) = temptarget(i);
          Oddcount = Oddcount + 1;
      end
  end
  target = [sort(Even) sort(Odd)]';
  clear Evencount Oddcount temptarget Even Odd
  minsep = min(abs(diff(target)));

 %end of ordering target entries
  sigma = norm(target - spectrum);
  sigma_0 =sigma;
  fprintf('\n  Starting with rho = %3.2f\n', rho)
  if exist('fileprint') == 1 fprintf('ties.out','\n  Starting with rho = %3.2f\n', rho),end

%Beginning of iteration:
  count = 0;
  countrho = 0;
  newstart = 1;  %to avoid recomputation of the spectrum of
     %toeplitz(tvec0)  and the associated Jacobian

  if rho > 0                        %We accept a solution from tiessub when
       subquit = min([10^(-jump1)*sigma,.1]);   %delta < subquit.
  else
       subquit = endquit;
  end

  fprintf('\n')
  if exist('fileprint') == 1 fprintf('ties.out','\n'),end
  while sigma >= 10^(-linquit)
      if rho >= .99999 break, end
      if badcount == badmax break,end
      target_rho = (1 - rho)*target+ rho*spectrum;
      tiessub
      if delta < subquit
           if exist('linmode') == 1 tvecold = tvec; end
           badcount = 0;
           sigma = norm(target - spectrum);
           if rho > 0
               subquit = 10^(-jump1)*sigma;
               fprintf('  ||spectrum - target|| = %3.2e  rho = %3.2f\n\n',sigma,rho)
               if exist('fileprint') == 1 fprintf('ties.out','  ||spectrum - target|| = %3.2e  rho = %3.2f\n\n',sigma,rho),end
            end
      end
      if  badcount == badmax
           badcount = 0;
           newstart = 1;
           spectrum = spectrum0;
           parity = parity0;
           rho = rho + deltarho;
           countrho = count;
           fprintf('  Increasing rho to %3.2f\n\n',rho)
           if exist('fileprint') == 1 fprintf('ties.out','  Increasing rho to %3.2f\n\n',rho),end
           subquit = 10^(-jump1)*sigma_0;
      end
      if badcount == badmax  badcount = 0;break,end
  end
  if exist('linmode') == 0
  if rho > 0  & rho < .99999
       badcount = 0;
       tvecsav = tvec;
       fprintf('  Setting rho = 0.00\n\n')
       if exist('fileprint') == 1 fprintf('ties.out','  Setting rho = 0\n\n'),end
   end
   if sigma >= endquit  & rho < .99999
       subquit = endquit;
       target_rho = target;
       rholast = rho; %save the value of rho in case terminal computation
       rho = 0;    %doesn't converge
       rhoswitch = 1;
       tiessub
   end
   end % if exist('linmode') == 0
if rho < .99999
if (delta  < subquit) | exist('linmode') == 1

  %Compute toepgen from normalized tvec. (The spectrum of toeplitz(tvec)
  %is target; the spectrum of toeplitz(toepgen) is target0.)

   sum1 = sum(target0);
   sum2 = norm(target0);
   toepgen = (sum2^2-sum1^2/n)^(1/2)*tvec;
   toepgen(1) = sum1/n;

   sorttarget = sort(target);

   if exist('screenprint') == 1
       fprintf('        target              tvec0       parity0        tvec       parity  \n\n')
       for i = 1:n
           fprintf('  %+15.14f   %+15.14f    %g  ',sorttarget(i),tvec0(i),parity0(i))
           fprintf('  %+15.14f    %g\n',tvec(i),parity(i))
       end
   end
   fprintf('\n')
   if exist('fileprint') == 1
       fprintf('ties.out','        target              tvec0       parity0        tvec       parity  \n\n')
       for i = 1:n
           fprintf('ties.out','  %+15.14f   %+15.14f    %g  ',sorttarget(i),tvec0(i),parity0(i))
           fprintf('ties.out','  %+15.14f    %g\n',tvec(i),parity(i))
       end
       fprintf('ties.out','\n')
   end
   clear sorttarget
   fprintf('ties.out','\n   %g  iterations with rho = %3.2f rhocount = %g\n',count,rholast, count - countrho)
 else % goes with "if delta < subquit
      fprintf(' \n Terminal computation did not converge\n')
      fprintf('  It was initiated with rho = %3.2f and I = %g\n', rholast, linquit)
      fprintf('\n')
      fprintf('ties.out',' \n Terminal computation did not converge\n')
      fprintf('ties.out','  It was initiated with rho = %3.2f and I = %g\n', rholast, linquit)
      fprintf('ties.out','\n')
 end  % goes with "if (delta < subquit) | exist('linmode')==1"

 else    % goes with "if rho  < .99999"
   fprintf('  rho = %3.2f. Too large!\n\n', rho)
   fprintf('\n')
   fprintf('ties.out','  rho = %3.2f. Too large!\n\n', rho)
   fprintf('ties.out','\n')
 end
 fprintf('  ||spectrum - target|| = %3.2e  rho =%3.2f', delta,rholast)
 fprintf('  count = %g   rhocount = %g\n\n', count, count - countrho)
 clear  A B P Q
 end % if exist('plotonly')==0

 if exist('plotonly')==1  plotpoly=1; end
 if exist('plotpoly')==1
    fprintf('\n  Preparing to graph results.\n')
    m=5*n;
    A = zeros(m,n);
       for rr = 1:m
           for ss = 1:n
              A(rr,ss) = (rr-1)*(ss-1);
           end
        end
    A = (pi/(m-1))*A;
    A = cos(A);
    y =2*A*tvec-tvec(1)*ones(m,1);
    x = 0:1/(m-1):1;
    plot(x,y,'b-')
    hold on
    for rr = 1:r
        plot([0,1],[E(rr,rr) E(rr,rr)],'g--')
    end
    for ss=1:s
        plot([0,1],[O(ss,ss) O(ss,ss)],'r--')
    end
    xlabel('Even eigenvalues in green, odd in red')
    hold off
 end
 if exist('plotonly')==1 clear plotonly plotpoly, end



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%Program tiessub.m
%tiessub.m this is a subroutine for tie.m

 if rhoswitch == 0

     delta = 25;
 else
     delta = sigma;
 end
 while  delta > subquit & badcount < badmax%beginning of major while loop
 delta_last = delta;

 if rhoswitch == 0 & newstart == 0 %This is the normal mode; however in switching from
     %\rho > 0 to terminal computation it is not necessary to recalculate sigma
 % If n is even go this way.
 A = zeros(r,r);
 B = zeros(s,s);
 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 = [(flipud(P0))',P0']';  %corresponding eigenvectors
     [Q0,O] = eig(A-B);       %odd eigenvectors
     Q = [-(flipud(Q0))',Q0']';   %corresponding eigenvectors
 else

 % If n is  odd go this way.

     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);  clear A %even eigenvalues
     Q0=P0;
     Q0([1],:)=[];
     P = [(flipud(Q0))',2^(1/2)*P0([1],:)',Q0']';
                    %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

     [P0,O] = eig(B);  clear B  %odd eigenvalues
     Q = [-(flipud(P0))',zeros(s,1),P0']';
                    %corresponding eigenvectors
 end

 [Even, order] = sort(diag(E));
 P=P(:,order);
 [Odd, order] = sort(diag(O));
 Q=Q(:,order); clear order

 %This finishes the computation of the eigenvectors.

 spectrum = [Even',Odd']'; %column vector with even eigenvalues in nondecreasing
                 % order, followed by odd eigenvalues in nondecreasing order.

 clear Even Odd  P0 Q0
 end  % goes with "if rhoswitch == 0 & newstart == 0"

 delta= norm(target_rho-spectrum)/norm(target_rho);

 if rhoswitch == 0
      if delta > delta_last  badcount = badcount + 1; end
 end
 if rhoswitch == 1 rhoswitch = 0;   end

 if rho == 0
     fprintf('  ||spectrum - target|| = %3.2e   ', delta)
     if exist('fileprint') == 1 fprintf('ties.out','  ||spectrum - target|| = %3.2e   ', delta),end
 else
     fprintf('  ||spectrum - target_rho|| = %3.2e   ', delta)
     if exist('fileprint') == 1 fprintf('ties.out','  ||spectrum - target_rho|| = %3.2e   ', delta),end
 end
 if badcount == badmax break,end

 if delta >= subquit
  if newstart == 0
     R = [P, Q];  clear P Q   % matrix of eigenvectors of $T$, each with norm=\sqrt2
     C = zeros(n,n);
     for i=1:n   C(i,1)=1; C(i,n) = R(1,i)*R(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) = 2*R(j:j+mm-1,[i])'*R(1:mm,[i])+R(j+mm,i)*R(mm+1,i);
                  end
         else
                  for i = 1:n
                      C(i,j) = 2*R(j:j+mm-1,[i])'*R(1:mm,[i]);
                 end
        end
     end
    end  % goes with "if newstart == 0"
  if newstart == 1
     C = C0;
     newstart = 0;
   end
     tvec= C\target_rho; clear C R
     count = count + 1;
     fprintf('count = %g \n', count)
     if exist('fileprint') == 1 fprintf('ties.out','count = %g  \n', count),end
 end

 end % of major "while" loop
 fprintf('\n\n')
 if exist('fileprint') == 1 fprintf('ties.out','\n\n'),end



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

% newties.m This program starts "ties.m using  Laurie's starting vector.
%(See D. P. Laurie, Initial values for the inverse Toeplitz eigenvalue
%problem, SIAM J. Sci. Comput., 22 (2001), 2239-2255.) To run it,
%enter target0 as a column vector. Then enter "newties."
%We highly recommend it; it's a winner!

t=toepinit(target0); % "toepinit.m is the next file
n=length(target0);
tvec0=t';
tvec0(1)=0;
s=0;
for r=1:n-1
   s=s+(n-r)*(tvec0(r+1))^2;
end
tvec0=tvec0/sqrt(2*s);
ties


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%Except for a minor adaptation, this function was created by Dirk P.Laurie.
function [t,M,Q]=toepinit(lam)

n=length(lam); lam=sort(lam(:))';
db=[0 filter(1,[1 1], lam)];
kc=[n+1:-2:1 rem(n+1,2)+2:2:n]; c=ifft(db(kc));
w=exp(pi*i*(0:n-1)/n);
ks=[n:-2:1 rem(n,2)+2:2:n]; s=w.*ifft(db(ks));
t=real(c+s); if nargout<2, return; end
[Q,D]=eig(toeplitz(t)); [d,k]=sort(diag(D));
Q=Q(:,k);
M=Q*diag(lam)*Q';

