%%%%generate random ... generically NOT hard case???
%%% TRS min x'Ax-2a'x  s.t. x'x <= s2
%%% TRS min (1/2)x'(2A)x-(2a)'x  s.t. x'x <= s2
%%% SDP: min   trace Abar*Xbar   s.t. trace Xbar <= 1+s2, Xbar_11 = 1
%%
clear all
n=10;
for ii=1:1,
%d=blkdiag(10*diag(rand(n,1)),-12,-10,-5*diag(rand(n-2,1)));
d=blkdiag(10*diag(rand(n,1)),-10.1,-10,-5*diag(rand(n-2,1)));
n=2*n;
%[q,~]=qr(randn(n));
%A=q*d*q';
%A=(A+A')/2;
A=d;
%%%Xbar=[x_0 x'
%%%      x X];
a=randn(n,1);
Abar=[0 -a'
      -a A];
s2=100*(rand+.1);   % tr radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%clear all
%%n=2;
%%A=[-1 0;0 1];
%%x=(sqrt(2)/2)*[1;1];
%%a=A*x;
%%s2=1;
n1=n+1;
In=speye(n);
AA=sparse([-In 2*A;2*A -4*a*a'/s2]);
BB=sparse([zeros(n) -In;-In zeros(n)]);
[va,da]=eig(2*A);
lb=-da(2,2);
ub=-da(1,1);
[lb ub]
[xv,lmb,iresult]=sptarn(AA,BB,lb,ub);
if isreal(lmb),
   nerror=norm(AA*xv-BB*xv*diag(lmb));
   if nerror>1e-10,
          fprintf('norm eigenvec error is %g \n',nerror)
   end       
   y1=xv(1:n,2);
   y2=xv(n+1:end,2);
   xs=-s2*xv(1:n,2)/((-2*a)'*xv(n+1:end,2));
   if abs(norm(xs)^2-s2)>1e-10,
   	fprintf('error in norm xs  %g \n',abs(norm(xs)^2-s2))
   end
   Lagrerror=norm((2*A+lmb(2)*In)*xs -2*a);
   if Lagrerror>1e-11,
   	fprintf('error in stationarity %g \n',norm((2*A+lmb(2)*In)*xs -2*a))
   end
   Vxs=null(xs');
   VAV= Vxs'*(2*A+lmb(2)*In)*Vxs;
   if min(eig(Vxs'*(2*A+lmb(2)*In)*Vxs)) <1e-12,
   	fprintf('2nd order cond not strict %g \n',...
	           min(eig(Vxs'*(2*A+lmb(2)*In)*Vxs)) )
   	keyboard
   end
fprintf('errors are \n')
[nerror abs(norm(xs)^2-s2) norm((2*A+lmb(2)*In)*xs-2*a) min(eig(Vxs'*(2*A+lmb(2)*In)*Vxs))] 
end   % end check lmb real

fprintf('problem %i  \n ',ii)
end   %  end of for loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%generate a loop for solving with increasing s2 to test Lemma 3.4 Martinez
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%for s2=s2:1:100*s2,
%%%cvx_begin  quiet sdp
%%%   cvx_precision best
%%%   variable Xbar(n1,n1) symmetric;
%%%   dual variable lambda;
%%%   minimize trace(Abar*Xbar);
%%%   subject to
%%%            Xbar(1,1)-1==0;
%%%            trace(Xbar)-1 - s2 <= 0: lambda;
%%%	    Xbar>=0;
%%%cvx_end
%%%fprintf(' ')
%%%xtrs=Xbar(2:end,1);
%%%[xv,~]=eig(Xbar);
%%%xtrse=xv(:,end);
%%%xtrse=xtrse/xtrse(1);
%%%xtrse=xtrse(2:end);
%%%pfeas=norm(xtrs)^2-s2;
%%%pfease=norm(xtrse)^2-s2;
%%%if max(pfeas,pfease) > 1e-11,
%%%	fprintf('error in pfeas value \n')
%%%	keyboard
%%%end
%%%pval=xtrs'*A*xtrs-2*a'*xtrs;
%%%pvale=xtrse'*A*xtrse-2*a'*xtrse;
%%%dfeas=norm( (A+lambda*eye(n))*xtrs -a)/(norm(a)+1);  %%relative stationarity 
%%%dfease=norm( (A+lambda*eye(n))*xtrse -a)/(norm(a)+1);  %%relative stationarity 
%%%if max(dfeas,dfease) > 1e-6,
%%%	fprintf('error in rel. dfeas value \n')
%%%	keyboard
%%%end
%%%fprintf('\n \n starting with new trs radius s2   %g  \n',s2)
%%%fprintf('primal feasibility ||xtrs||^2-s2 %g  \n',pfeas)
%%%fprintf('rel. dual feasibility/stationarity %g \n',dfeas)
%%%fprintf('pval(x),pval(Xbar)  %g %g  \n',pval,trace(Abar*Xbar))
%%%if abs(pval-trace(Abar*Xbar))/(abs(pval)+1)>1e-8,
%%%	fprintf('error in pval values \n')
%%%	keyboard
%%%end
%%%mineig=min(eig(A+lambda*eye(n)));
%%%fprintf('min eig Hessian Lagr %g \n',mineig)
%%%if mineig<0,
%%%	fprintf('error in mineig %g \n',mineig)
%%%	keyboard
%%%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%for shift=.1:.01:.99,
%%%%%%%now add a linear constraint to cut off global min. xtrs
%%%%%%  let b=(1-shift)*xtrs and alpha= norm(b)^2
%%%%%% equation is 0 <= alpha-b'x  note that alpha =shift^2*delta 
%%%%%%                 so xtrs is infeas.
%%%%%% Add the second order cone constraint
%%%%%%  norm(alpha*x-X*b)<=alpha-b'x
%%%b=shift*xtrs;
%%%alpha=norm(b)^2;
%%%cvx_begin quiet sdp
%%%   cvx_precision best
%%%   %variable Xbar(n1,n1) symmetric;
%%%   variable X(n,n) symmetric;
%%%   variable x(n);
%%%   dual variable muu;
%%%   minimize trace(A*X)-2*(a'*x);
%%%   subject to
%%%            %Xbar(1,1)-1==0;
%%%            %trace(Xbar)-1 - s2 <= 0: lambda;
%%%            trace(X) - s2 <= 0: muu;
%%%	    norm(alpha*x-X*b)-(alpha-b'*x) <= 0;
%%%	    %Xbar>=0;
%%%	    [1 x'
%%%	      x X ] >=0;
%%%cvx_end
%%%fprintf(' ')
%%%pfeas=norm(x)^2-s2;
%%%pfease=norm(x)^2-s2;
%%%if max(pfeas) > 1e-11,
%%%	fprintf('error in pfeas value \n')
%%%	keyboard
%%%end
%%%pval=x'*A*x-2*a'*x;
%%%fprintf('\n \n starting with new shift   %g  \n',shift)
%%%fprintf('primal feasibility ||x||^2-s2 %g  \n',pfeas)
%%%mineig=min(eig(A+lambda*eye(n)));
%%%fprintf('min eig Hessian Lagr %g \n',mineig)
%%%if mineig<0,
%%%	fprintf('error in mineig %g \n',mineig)
%%%	keyboard
%%%end
%%%lam=sort(eig(A+muu*eye(n)),'ascend');
%%%fprintf('lam1 muu lam2 for Hessian Lagr %g  %g  %g \n',lam(1),muu,lam(2))
%%%linineq=alpha-b'*x;
%%%fprintf('value of linear inequality - test *active*  %g  \n',linineq)
%%%if linineq>1e-3,
%%%	fprintf('lin constr NOT active \n')
%%%	keyboard
%%%end
%%%end   %   end for 'shift'
%%%end   %   end for s2 increases
