%% function [cvx_optval,X,Y,W] = Diamond(J,dim)
% computes the diamond norm of the Hermiticity-preserving map J
% with input/output dimensions optionally provided in dim=[dimA,dimB]

function [cvx_optval,X,Y,W] = Diamond(J,dim)
    
J = (J+J')/2; % make sure that J is Hermitian, avoid numerical issues

if nargin<2
    d = max(size(J));
    dim = [1,1]*sqrt(d);
else
    d = dim(1)*dim(2);
end
 
cvx_begin sdp quiet

   variable X(d,d) hermitian semidefinite
   variable Y(d,d) hermitian semidefinite
   variable mu nonnegative
   dual variable W
   
   PartialTrace(X+Y,2,dim) <= mu*eye(dim(1));
   W: X-Y == J;
   
   minimize mu
   
cvx_end
   
end