function [prob,C,G,X,Y,Z,lam] = prob_mio(rho, m, tol)
% prob_mio(rho,m,tol) computes the maximal probability of distillation of an m-dimensional maximally coherent state under DIO with tolerance tol
% rho: dxd density matrix or normalised state vector
% m: target state dimension
% tol: error tolerance 0<=tol<=1

d=size(rho,2);
if d==1
    rho=rho*rho';
    d=size(rho,1);
end
rho=(rho+rho')/2;

cvx_begin sdp quiet

    variable C(d,d) hermitian
    variable G(d,d) hermitian
    dual variable X
    dual variable Y
    dual variable Z
    dual variable lam
    
    maximize ( trace(G*rho) )

    0 <= C
    Y : C <= G
    X : G <= eye(d)

    Z : diag(diag(G))==m*diag(diag(C))

    lam : trace(C*rho) >= (1-tol)*trace(G*rho)

cvx_end

prob = cvx_optval;

end