function [U,V] = symX(U,V,RealFlag)
% Enforces symmetry of X=V/U with [U;V] orthogonal.

%   Copyright 2018 The MathWorks, Inc.

% Refs:
%  * QZ: Laub and Gahinet, TAC, 42 (1997), pp. 1303-1307.
%  * Polar: Higham, Matrix Nearness Problems and Applications
% Note:
%  * While polar approach causes smaller rotation of [U;V], it can lead
%    to a much larger asymmetry in V/U.

if RealFlag
   [DU,DV,Qr,~] = qz(U,V,'real');
   DU = diag(DU);  DV = diag(DV);
   DS = sqrt(DU.^2+DV.^2);
   U = Qr'*diag(DU./DS);
   V = Qr'*diag(DV./DS);
   % Polar:
   %    [P,~,Q] = svd(complex(U,V));
   %    PQ = P*Q';
   %    U = real(PQ);  V = imag(PQ);
else
   [DU,DV,Qr,~] = qz(U,V);
   DU = diag(DU);  DV = diag(DV);
   % Enforce |DU|^2+|DV|^2=1 and DV/DU real
   theta = atan2(real(DU.*conj(DV)),(abs(DU).^2-abs(DV).^2)/2)/2;
   c = cos(theta);  s = sin(theta);
   tau = c.*DU+s.*DV;  tau = tau./abs(tau);
   U = Qr'*diag(c.*tau);
   V = Qr'*diag(s.*tau);
end

% [~,sa] = ltipack.angleS([U0;V0],[U;V]);
% X = V/U;
% QZ = [norm(U0'*V0-V0'*U0) norm(U'*V-V'*U) norm(X-X',1)/norm(X,1) max(sa)]
   

