function capacity = symcap(H,P,its,tol)
% Symmmetric Capcity
%
%   Usage:  cap=symcap(H,P)
%           cap=symcap(H,P,its)
%           cap=symcap(H,P,its,tol)
%   Inputs: 
%           H : 3-dim channel matrices
%           P : power constraint
%           its : dpc steps
%           tol : tolerance for the ellipsoid algoirithm
%   Output: 
%           cap : symmetric rate for an individual user
%
% Authors: Juyul Lee and Nihar Jindal
% Date: 05/31/2006
%
% Note: This routine is based on the paper:
%       J. Lee and N. Jindal
%       "Symmetric capacity of a downlink MIMO channel,"
%       IEEE ISIT 2006.

if nargin < 4
    tol = 1e-4;
end;

if nargin < 3
    its = 25;
end;

[M,N,K]=size(H);

mu = (1/K)*ones(K,1); % init
E = (1-1/K)*eye(K-1);
subgrad=zeros(K-1,1);

while(1)
    [rate_y,Q]=f_wtdratesumcap2(mu,H,10*log10(P),its);
    
    % check
    sort_rate = sort(rate_y);
    mean_rate = mean(rate_y);
    if((mean_rate-sort_rate(1)<tol) &&(sort_rate(K)-mean_rate<tol))
        break;
    end;
    
    if(max(eig(E)) < 0.001*tol)
        break;
    end;
    
    % subgradient
    for k=1:K-1
        subgrad(k) = rate_y(k)-rate_y(K);
    end;
    subgrad_n = subgrad / sqrt(subgrad'*E*subgrad);

    % Updates
    if(K==2)    % special case of the ellipsoid algorithm
        mu(1:K-1) = mu(1:K-1) - (1/K)*E*sign(subgrad_n);
        E=0.5*E;
    else
        mu(1:K-1) = mu(1:K-1) - (1/(K))*E*subgrad_n;
        E = ((K-1)^2/((K-1)^2-1))*(E-(2/(K))*E*subgrad_n*subgrad_n'*E);
    end;
    mu(K) = max(1 - sum(mu(1:K-1)),0);
    
    if(length(find(eig(E)<0))>0)
        break;
        %error('Negative length of axis due to finite representation...');
    end;
end;

capacity = sum(rate_y)/K;