%% User:   can change: dimension, sparsity for random sparse Q 
%%       global variables must be used
clear all
profile on
global Xu Q q2kron indsu ;
itercnt=[];
dimcnt=[];
densitycnt=[];
normZXcnt=[];
nnzQcnt=[];
mineigcnt=[];
cpucnt=[];
optvalcnt=[];
crossitercnt=[];
crossindex=1;
%for n=365:365,
 for n=35:35,
%n=35;for crossindex=1:2:12,
%n=45;for ii=1:6,
  %density=1/(2*(ii^3));
  density=1/(n);         % density decreases with n
  tforrandQ=clock;
  Q=5*sprandsym(n,density);
  e=ones(n,1);
  r=9*randn(n,1);
  r2=[0;r(1:n-1)];
 %Q=spdiags([r 0*e r2], -1:1,Q);  % r2 used due to 'bug' in spdiags???
                                  % i.e. r2(1) is skipped  by spdiags
  if norm(full(Q)-full(Q)','fro')>.000000001,
   'Q symm error'
     keyboard
   end
  Q=spdiags([e 0*e e], -1:1,Q);
  Q=round(Q);
  randQtime=etime(clock,tforrandQ);
  disp(['time for generating random Q:  ',num2str(randQtime)]);
  nnzQ=nnz(Q);
  nnzQcnt=[nnzQcnt nnzQ];
  disp(['number nonzeros in  Q:  ',num2str(nnzQ)]);
  [m,n]=size(Q);
%%%%%%%%%%%%%%%%%%call to solve max-cut using sdp relax%%%%%%%%%%%%%%
  toler=10^(-14);
  crosstoler=10^(-crossindex);
  disp(['solve random Q dim, crosstoler :  ',num2str(n),...
          '  ',num2str(crosstoler)]);
  tcpusolve=cputime;
  [X,y,Z,mu,itersinfo,meritfn]=mcsdp(toler,crosstoler);
  tcpusolve=cputime-tcpusolve;
%%%%%%%%%%%%%%%%%%%%%% end of call to MCSDP
  mineig=min(min(eig(X)),min(eig(Z)));
  mineigcnt=[mineigcnt mineig];
  itercnt=[itercnt itersinfo(1)];
  crossitercnt=[crossitercnt itersinfo(2)];
  dimcnt=[dimcnt n];
  normZXcnt=[normZXcnt norm(Z*X,'fro')];
  densitycnt=[densitycnt density];
  cpucnt=[cpucnt tcpusolve];
  optvalcnt=[optvalcnt trace(Q*X)];
  save sdpdata itercnt dimcnt densitycnt nnzQcnt mineigcnt cpucnt optvalcnt normZXcnt n crossitercnt
profile report
end