function [User, max_diff] = PowerMinAltMthd(User,Cells,Chn,BS,Sim,intf_2nd)

%%%% This function minimizes the total transmit power under fixed Rx BF
%%%% subject to
%%%% per-BS power constraints and per-user rate constraints
ant_gain = 1;
max_diff=0;

Num_Cell = Sim.Num_Cell;
Total_BS = Sim.Num_BS;
L_Macro = Sim.Num_Macro;
L_Pico = Sim.Num_Pico;
Num_TxAnt_Macro = Sim.Num_Ant_Macro;
Num_TxAnt_Pico = Sim.Num_Ant_Pico;
Num_Rx_Ant = Sim.Num_Rx_Ant;
Q = Num_Rx_Ant;


fprintf('$$$$$ Power Minimization Procedure $$$$$\n');

[User, Total_Ant,Per_BS_Coop,Per_BS_ServAnt] = Compute_Coop_Link(Num_Cell,User,Cells,Num_TxAnt_Macro,Num_TxAnt_Pico,L_Macro,L_Pico);

N = Total_Ant;
previous_sumpower = 0;
iter = 1;
while iter <= 5
    
    %%%%%%%%%%%%%%%%%%%%
    %compute MMSE Rx BF
    %%%%%%%%%%%%%%%%%%%%
    for l = 1:Num_Cell
        for k_index = 1:length(Cells(l).Scheduled_User)
            k = Cells(l).Scheduled_User(k_index);
            H = Get_Chn(Chn, User(l,k).ServingCluster, l,k, Sim);
            
            Kz_DL = eye(Q)* User(l,k).eff_rx_noise + intf_2nd;
            %total received signal covariance matrix for user k
            
            for m = 1:Num_Cell
                for kk_index = 1:length(Cells(m).Scheduled_User)
                    kk = Cells(m).Scheduled_User(kk_index);
                    beam = User(m,kk).beam_tx;
                    Hkk = Get_Chn(Chn, User(m,kk).ServingCluster, l,k, Sim);
                    Kz_DL = Kz_DL + ant_gain*Hkk'*(beam*beam')*Hkk;
                end
            end
            
            beam_rx = Kz_DL\(H'*User(l,k).beam_tx)*sqrt(ant_gain);
            User(l,k).beam_rx = beam_rx;
        end
    end
    
    %%%%%%%%%%%%%%%%%
    %Optimize Tx BF
    %%%%%%%%%%%%%%%%%
    
    %fprintf( 'call cvx solver \n');
    fprintf('$$$$ First try sdpt3 $$$$\n');
    cvx_begin quiet
    cvx_solver sdpt3;
    variable V(N,1) complex; %Tx BF
    
    minimize( norm(V) );
    subject to
    
    %per BS power constraint
    for m = 1:Num_Cell
        for l = 1:Total_BS
            
            if ~isempty(Per_BS_Coop(m,l).Coop_Matrix) && ~isempty(Per_BS_ServAnt(m,l).index)
                norm(Per_BS_Coop(m,l).Coop_Matrix * V(Per_BS_ServAnt(m,l).index)) <= sqrt(BS(m,l).bf_power_constraint); %per BS power constraint
            end
            
        end
    end
    
    %per user rate constraint
    ant_idx = 1;
    for l = 1:Num_Cell
        for ik = 1:length(Cells(l).Scheduled_User)
            k = Cells(l).Scheduled_User(ik);
            
            gamma_k = (2^(User(l,k).inst_rate) - 1)*Sim.snr_gap; %SINR target
            
            if gamma_k <= 0
                fprintf('!!!!!!! User SINR is NOT CORRECT !!!!!!!\n');
                return;
            end
            
            h_k_blk_diag = [];
            for m = 1:Num_Cell
                for ikk = 1:length(Cells(m).Scheduled_User)
                    kk = Cells(m).Scheduled_User(ikk);
                    
                    h_jk = (User(l,k).beam_rx)' * Get_Chn(Chn, User(m,kk).ServingCluster, l, k, Sim)';   %row vector
                    h_k_blk_diag = blkdiag(h_k_blk_diag, h_jk);
                    
                    if (m == l) && (k == kk)
                        h_kk = h_jk;
                    end
                    
                end
            end
            
            ant_range = ant_idx : (ant_idx+User(l,k).ServingAnt-1);
            
            h_kk*V(ant_range) == real(h_kk*V(ant_range));
            
            norm( [h_k_blk_diag*V; sqrt(User(l,k).eff_rx_noise)*norm(User(l,k).beam_rx)] ) <= 1/2 * sqrt(1 + 1/gamma_k) * ( h_kk*V(ant_range) + (h_kk*V(ant_range))');
            
            ant_idx = ant_idx + User(l,k).ServingAnt;
        end
    end
    cvx_end
    cvx_status
    
    if ~(strcmp(cvx_status, 'Solved') || strcmp(cvx_status, 'Inaccurate/Solved'))
        fprintf('$$$$ Second try sedumi $$$$\n');
        cvx_begin quiet
        cvx_solver sedumi;
        variable V(N,1) complex; %Tx BF
        
        minimize( norm(V) );
        subject to
        
        %per BS power constraint
        for m = 1:Num_Cell
            for l = 1:Total_BS
                
                if ~isempty(Per_BS_Coop(m,l).Coop_Matrix) && ~isempty(Per_BS_ServAnt(m,l).index)
                    norm(Per_BS_Coop(m,l).Coop_Matrix * V(Per_BS_ServAnt(m,l).index)) <= sqrt(BS(m,l).bf_power_constraint);  %per BS power constraint
                end
                
            end
        end
        
        %per user rate constraint
        ant_idx = 1;
        for l = 1:Num_Cell
            for ik = 1:length(Cells(l).Scheduled_User)
                k = Cells(l).Scheduled_User(ik);
                
                gamma_k = (2^(User(l,k).inst_rate) - 1)*Sim.snr_gap; %SINR target
                
                if gamma_k <= 0
                    fprintf('!!!!!!! User SINR is NOT CORRECT !!!!!!!\n');
                    return;
                end
                
                h_k_blk_diag = [];
                for m = 1:Num_Cell
                    for ikk = 1:length(Cells(m).Scheduled_User)
                        kk = Cells(m).Scheduled_User(ikk);
                        
                        h_jk = (User(l,k).beam_rx)' * Get_Chn(Chn, User(m,kk).ServingCluster, l, k, Sim)';   %row vector
                        h_k_blk_diag = blkdiag(h_k_blk_diag, h_jk);
                        
                        if (m == l) && (k == kk)
                            h_kk = h_jk;
                        end
                        
                    end
                end
                
                ant_range = ant_idx : (ant_idx+User(l,k).ServingAnt-1);
                
                h_kk*V(ant_range) == real(h_kk*V(ant_range));
                
                norm( [h_k_blk_diag*V; sqrt(User(l,k).eff_rx_noise)*norm(User(l,k).beam_rx)] ) <= 1/2 * sqrt(1 + 1/gamma_k) * ( h_kk*V(ant_range) + (h_kk*V(ant_range))');
                
                ant_idx = ant_idx + User(l,k).ServingAnt;
            end
        end
        cvx_end
        cvx_status
        
        if ~(strcmp(cvx_status, 'Solved') || strcmp(cvx_status, 'Inaccurate/Solved'))
            fprintf('$$$$ Both sdpt3 and sedumi do not work $$$$\n');
            break;
        end
    end
    
    % if strcmp(cvx_status, 'Infeasible')
    %     cvx_status
    % end
    
    %assign cvx solution to beam_tx and find the max beam_tx difference for convergence check
    max_diff = 0;
    ant_head = 1;
    for l = 1:Num_Cell
        for ik = 1:length(Cells(l).Scheduled_User)
            k = Cells(l).Scheduled_User(ik);
            
            %previous_beam = User(l,k).beam_tx;
            ant_tail = ant_head + User(l,k).ServingAnt -1;
            
            User(l,k).beam_tx = V(ant_head:ant_tail);
            %diff = norm(User(l,k).beam_tx - previous_beam)/norm(previous_beam);
            
            %if diff > max_diff
            %    max_diff = diff;
            %end
            
            ant_head = ant_tail+1;
            
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%
    %check convergence
    %%%%%%%%%%%%%%%%%%%%
    
    current_sumpower = norm(V)^2;
    power_diff = (current_sumpower - previous_sumpower)/current_sumpower;
    fprintf('$$$$$ iter = %d, sum power diff = %f $$$$$\n',iter, power_diff);
    if abs(power_diff) < 1e-2
        break;
    else
        iter = iter + 1;
        previous_sumpower = current_sumpower;
    end
end
