function X = parrott(B,C,D,GAM)
% Solves Parrott completion problem.
%
%   Given matrices B,C,D and 
%
%     G >=  GMIN := MAX ( || [B,D] || , || [C;D] || )
%
%   computes a matrix X such that
%
%        || X   B ||
%        ||       ||   <=  G
%        || C   D ||2

%   Copyright 2018 The MathWorks, Inc.
p1 = size(B,1);
m1 = size(C,2);
[p2,m2] = size(D);
if GAM==Inf || GAM==0
   X = zeros(p1,m1);
else
   % Note: Use SVD to handle the case ||D||=GAM
   %   * If || D*v || = GAM, then B*v=0 from || [B;D] || <= GAM
   %   * if || u'*D || = GAM, then u'*C=0 from || [C,D] || <= GAM
   % We can therefore ignore directions where [GAM*I D;D' GAM*I] is singular.
   nsv = min(p2,m2);
   [u,sD,v] = svd(D);
   s = diag(sD(1:nsv,1:nsv));
   NonSing = (s<(1-1e3*eps)*GAM);
   if all(NonSing)
      % [GAM*I D;D' GAM*I] invertible
      X = [zeros(p1,p2),B] * ([GAM*eye(p2) D;D' GAM*eye(m2)] \ [C;zeros(m2,m1)]);
   else
      ix = find(NonSing.');
      iu = [ix,nsv+1:p2];  nu = numel(iu);
      iv = [ix,nsv+1:m2];  nv = numel(iv);
      aux = sD(iu,iv);
      X = [zeros(p1,nu) , B*v(:,iv)] * ...
         ([GAM*eye(nu) aux;aux' GAM*eye(nv)] \ [u(:,iu)'*C ; zeros(nv,m1)]);
   end
end