%
%  This function obtains a unitary matrix Q such that: d=diag(Q'*diag(lmd)*Q). 
%  In other words, it gives a way to generate a matrix with given eigenvalues and diagonal elements.
%
%  By Daniel Perez Palomar (last revision: May 10, 2004).
%  Feel free to distribute this file as it is (without including any modifications).
%  Please, email any improvement or error to Daniel.P.Palomar@ieee.org or danielp@princeton.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Details of the function:
% ------------------------
%
% It uses the method in [Viswanath and Anantharam IT99] and also the method
% in [Marshall79] for verification purposes. See [Palomar-PhD2003] for details.
%
% INPUTS: Vectors of eigenvalues and diagonal elements, lmd and d, in
% decreasing order (of course, lmd must majorize d, otherwise there is no
% such unitary matrix Q).
%
% OUTPUT: Unitary matrix Q such that: d=diag(Q'*diag(lmd)*Q).
%
function Q_acc = Q_lmd_diag(lmd,d)

EPS=1e-12;  %small numerical error allowed in the comparisons
N=length(lmd);
if any(flipud(sort(lmd))~=lmd), error('lmd is not in decreasing order !!'); end
if any(flipud(sort(d))~=d), error('d is not in decreasing order !!'); end
if abs(sum(lmd)-sum(d))>EPS, error('sum(lmd)!=sum(d) !!'); end
for n=1:N, if sum(lmd(1:n))<sum(d(1:n))-EPS, error('lmd does not majorize d !!'); end; end % if sum(lmd(1:n)) < sum(d(1:n))

mask=ones(N); for i1=1:N, for i2=1:i1-1, mask(i1,i2)=-1; end; end 
T_acc=eye(N);
Q_acc=eye(N);
n=N;
while 1,
   %
   %Find both indexes
   %
   while n>=1 && abs(lmd(n)-d(n))<EPS, n=n-1; end   %if n>1 & lmd(n)==d(n)
   if n<=1, break; end; %finish method
   k=n-1;
   while lmd(k)-d(k)<=EPS, k=k-1; if k<1, error('k==0 in Q_lmd_diag !!'); end; end   %if lmd(k)<=d(k)
   %
   % Find alpha
   %
     %Find alpha so that: (1-alpha)*lmd(k) + alpha*lmd(n) = d(n)
     alpha = (d(n)-lmd(k))/(lmd(n)-lmd(k));
     %check whether this alpha is too large
     if alpha*lmd(k) + (1-alpha)*lmd(n) < d(k),
   	  alpha = (d(k)-lmd(n))/(lmd(k)-lmd(n));
     end
     %Using the algorithm in [Marshall79]
     alpha_bis=1-min(d(n)-lmd(n),lmd(k)-d(k))/(lmd(k)-lmd(n));  
     if abs(alpha_bis-alpha)>1e-12, error('Algorithms in [Marshall79] and [Viswanath99] differ !!'); end
   %  
   % T-transform for the transformaton of the vector
   %
   T=eye(N);
   T(k,k)=alpha; T(n,n)=alpha;
   T(k,n)=1-alpha; T(n,k)=1-alpha;
   lmd=T*lmd;
   T_acc=T*T_acc;
   %unitary matrix for the transformation of the matrix with diagonal elements
   Q_acc=Q_acc*(sqrt(T).*mask);
end