%%%%%%%% IntroExample %%%%%%%%%%%%
%
%Example from p.2 of the manuscript
%
%Sremac, Woerdeman, and Wolkowicz
%2020

n=5;
tn = n*(n+1)/2;     %triangular number

%%%%%%%%%%%%%%%%%%%%%for ssMat/ssvec
n2=n^2;
svind = find(triu(ones(n),1)); %linear indices of the strictly upper triangle
sMind = zeros(n2,1); %matrix indices with respect to linear indices of the diagonal and the strictly upper triangle
sMind(svind) = 1:tn-n; %strictly upper triangle
sMind(1:n+1:n2) = tn-n+1:tn;%the diagonal
%calculate the linear indices of the strictly lower triangle
indM = zeros(n); %matrix indices
indM(:) = 1:n2;
indM = indM';
lsvind = indM(svind); %strictly lower triangular
sMind(lsvind) = 1:tn-n; %strictly lower triangular

%%%%%%%%%%%%%%%%%%%%% make problem
Rot = eye(n);       %Orthogonal matrix, if any rotation/reflection is desired, none by default
SDlevent = 4;       %singularity degree desired
tracebound = 1;     %bound for trace to get bounded problem
[A,b,Xorig] = MakeProblemTun(n,tracebound,svind,SDlevent,Rot);   %construct A,b,Xorig
fprintf('\n Generating and solving the example at the beginning of the manuscript.\n')
fprintf('dim = %i,  sing. deg = %i (can be changed in file) \n',n,SDlevent)

%%%%%%%%%%%%%%%%%%%%% solve with cvx
cvx_clear
try
    cvx_solver mosek;   %using cvx with mosek, try with mosek first
catch
    warning('Problem using solver mosek.  Changing to SDPT3.');
    cvx_solver sdpt3;   
end
cvx_save_prefs
cvx_begin sdp
variable X(n,n) symmetric
minimize 0
subject to
A*ssvec(X,svind) == b;
X >= 0;
cvx_end

%%%%%%%%%%%%%%%%%%%%% errors
ferror = norm(X-Xorig,'fro');   %forward error
berror = norm(A*ssvec(X,svind)-b) + max(0,-min(eig(X))); %backward error
fprintf('forward error %g;  backward error %g\n',ferror,berror)
