function [x_opt, Q_opt, fval_opt, fgap_opt] = mutual_mimo_opt(H0, Rt, M, N, SNR_db)%----------------------------------------------------------------------% Optimization function to find the ergodic capacity of a MIMO channel given% the channel mean and transmit correlation matrices. The program returns% the optimum variable values (which can be the eigenvalues of the transmit% covariance matrix or the vectorized version of this matrix), the optimum% transmit covariance matrix itself, the mutual information value and the% function gap to optimality.%% Author: Mai Vu% Date: Dec 04, 2003% Modified: Mar 02, 2004; Apr 19, 2004.%----------------------------------------------------------------------%----------------------------------------------------------------------% INPUTS:%%   N: Number of transmit antennas%   M: Number of receive antennas%   H0: channel mean (MxN)%   Rt: transmit correlation%   SNR_db: SNR in dB%% THE PROBLEM%%   max f = E[logdet(I + SQ)]%   s.t  trace(Q) = 1%%   where  S = H'*H%          H = H0 + Hw*Rt_sqrt%% OUTPUTS:%   %   x_opt: In symmetric Q case, the vector contains re-order elements of Q%   Q_opt: Optimum input covariance matrix for each SNR input%   fval_opt: Mutual information value at each SNR input%   fgap_opt: Gap to optimality of the MI value at each SNR input%----------------------------------------------------------------------  % optimization parametersNSAMPLE = 10000;MAXITER = 10;MAXBARR = 10;MAXLINES = 50;epsilon = 5*10^-6;mu = 100;barr_t = 100;% signal to noise ratiogamma = 10.^(SNR_db/10);I = eye(N,N);if (H0 == 0)  % special case correlation only  [U Rt] = eig(Rt);  diag_flag = 1;elseif (Rt/trace(Rt)*N == I)  % special case mean only  [V H0 U] = svd(H0);  diag_flag = 1;else  diag_flag = 0;endRt_sqrt = sqrtm(Rt);% initial transmit covariance matrixS0 = H0'*H0 + M*Rt;Q = S0/trace(S0);lambda = real(eig(Q));% indices for calculating the gradient/hessian and reorganize the unknowns% into a matrixind1 = [];ind2 = [];for n=1:N  ind1 = [ind1, n:1:N];  ind2 = [ind2, 1:1:N-n+1];endind3 = [];for n=1:N  if n>1    for k=1:n-1      ind3 = [ind3 k+(n-k)*N-(n-k-1)*(n-k)/2 + N*(N-1)/2];    end  end  for k=n:N    ind3 = [ind3 n+(k-n)*N-(k-n-1)*(k-n)/2];  endend% number of unknownsNr = (N+1)*N/2;Nx = N^2;% initial unknown valuesQ12 = diag(Q(ind1,ind2));x = [real(Q12(1:N))/2; real(Q12(N+1:Nr)); imag(Q12(N+1:Nr))];% initial line search counter arrayls_cnt = [];% maximum number of barrier outter stepsbarr_cnt = ceil(log(1/epsilon/barr_t)/log(mu)) + 1;% solve for diagonal Qif (diag_flag)    for k=1:length(gamma)    disp(['  SNR = ', num2str(SNR_db(k)), 'dB']);    [lambda_new, fval2, fgap2, cnt2(k)] = ...        mutual_mimo_newton_diag(M, N, H0, Rt_sqrt, gamma(k), lambda, ...                                NSAMPLE, MAXITER, MAXLINES, epsilon);    x_opt(:,k) = lambda_new(:,cnt2(k));    Q_opt(:,:,k) = U*diag(x_opt(:,k))*U';    fval_opt(k) = fval2(cnt2(k));    fgap_opt(k) = fgap2(cnt2(k));  endelse% solve for symmetric Q  for k=1:length(gamma)    disp(['  SNR = ', num2str(SNR_db(k)), 'dB']);    [x_new, fval, fgap, cnt_temp, ls_cnt_new] = ...        mutual_mimo_newton_symm_choose(M, N, H0, Rt_sqrt, gamma(k), Q, NSAMPLE, ...                                       MAXITER, MAXLINES, epsilon, ind1, ind2, ...                                       ind3, barr_t, mu, MAXBARR);    x_opt(:,k) = x_new(:,cnt_temp(length(cnt_temp)));    Q_opt(:,:,k) = form_Q_symm(x_opt(:,k), N, Nr, Nx, ind3);    fval_opt(k) = fval(length(fval));    fgap_opt(k) = fgap(length(fgap));  endend